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MODEL-BASED  ARRAY  PROCESSING 


1.  INTRODUCTION 


Sonar  systems  are  being  driven  to  use  ever  lower  frequencies  to  compensate  for  the  quieting  of  acoustic 
emissions  from  targets  of  interest.  Because  these  low-frequency  sounds  are  more  difficult  to  suppress  at 
the  source  and  thus  can  propagate  further  through  the  underwater  medium,  reasonable  detection  ranges 
can  be  maintained. 

For  a  given  sonar  receiver  and  processing  chain,  the  beamwidth  of  the  system  increases  as  the  frequency 
decreases.  As  the  beams  become  broader,  there  is  a  greater  chance  that  more  than  one  acoustic  source  will 
be  present  in  the  same  beam.  This  report  introduces  beamforming  methods  that  are  aimed  at  detecting 
more  than  one  source  within  a  single  beam.  The  techniques  are  based  on  planewave  models  and  use  fitting 
procedures  to  find  optimal  sets  of  planewaves,  in  a  minimal  squared-error  sense,  to  match  measured  time 
series  data. 

Each  of  the  following  sections  first  presents  the  theoretical  development  and  then  discusses  the  results 
with  simulated  data,  comparing  these  results  with  conventional  beamforming  techniques. 

In  section  2,  a  single  or  one-planewave  model  is  fitted  to  data  observed  by  a  line  array  of  sensors  that 
can  be  arbitrarily  spaced;  also  discussed  is  the  special  case  of  a  sparse  equispaced  array. 

In  section  3,  a  two-planewave  model  is  fitted  to  linear  array  data  (arbitrarily  spaced).  Section  4  applies 
a  similar  technique  to  the  situation  where  the  two  planewaves  are  related  by  a  single  change  in  amplitude 
and  a  time  delay.  This  case  represents  the  arrival  of  two  planewaves  from  a  single  source  via  two  different 
propagation  paths  (multipath). 

Finally,  in  section  5,  a  model  is  fitted  of  a  planewave  arriving  at  an  arbitrarily  spaced  two-dimensional 
array  from  a  moving  source.  In  this  case,  the  fitting  procedure  yields  two-dimensional  start  and  end 
positions  of  the  source,  and  indicates  the  strengths  of  its  frequency  components. 

Model-based  processing  can  be  used  as  a  preprocessor  prior  to  such  additional  beamforming  as  con¬ 
ventional,  adaptive,  Fourier  Integral  Method  (FIM),  etc.  If  a  double  source  is  found,  the  stronger  one  can 
be  coherently  removed,  using  the  modeled  coefficients  yielded  by  the  fitting  procedure,  to  better  identify 
the  weaker  source(s).  Removed  sources  can  be  added  to  subsequent  sonar  displays  to  prevent  loss  of 
information  to  operators. 
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2.  ONE-PLANE  WAVE  FIT  FOR 
UNEQUALLY  SPACED  LINE  ARRAYS 


2.1  THEORY 

2.1.1  Arbitrarily  Spaced  Line  Array 

A  continuous  pressure  field  p(t,  x )  is  assumed  at  time  t  and  location  x.  Time  samples  are  taken  at 
t  =  nA  for  n  =  0:1V— 1,  and  spatial  samples  are  taken  at  locations  x  =  x(m)  for  m  =  0:M  —  1.  The 
colon  symbol  J:K  is  used  to  denote  {J,J+ 1, . . . , K}  with  J  <  K.  The  total  discrete  data  available  are 
p(n, m)  =  p(nA, x(m))  for  n  =  0 : IV - 1,  m  =  0:M - 1. 

If  a  single  planewave  arrives  from  angle  Ui  =  sin  8\  comprised  of  frequencies  {/i  (&)}£* ,  the  observed 
complex  pressure  field  at  ( t ,  x)  is  modeled  as 

k„ 

Pi  (t,  x)  =  ^  ai  (k)  exp  [i  2ir  fi  ( k)(t  -  f  ui )] ,  (1) 

k=ka 

where  the  hypothesized  amplitudes  {aj(A:)}*‘  are  complex.  For  given  frequencies  {/i(A:)}*‘  and  arrival 
angle  iq,  amplitudes  should  be  chosen  so  that  the  total  weighted  fitting  error  e  is  minimized. 

Here,  error  is  defined  as 


N- 1  M- 1 


data 


fit 


e=  J!  wt  (n) w*  (m)  |p(nA>  x(m))  -  pi  (nA,  x(m))  | 

n=0  m=0 

=  ^ wt(n)wx{m )  [p(n, m)  —  ^ai(A:)  expji  a\{k)  [n  —  /?(m)ui]  j 


n,m 


(2) 


where  the  temporal  and  spatial  weights,  {u>t (n)}  and  {wx(m)},  are  real  and  positive,  and  the  known 
dimensionless  parameters  are 


Qi(/c)  =  27r/i(fc)A,  0(m)  = 


To  minimize  e,  Kay’s  partial  derivative  procedure1  is  used;  thus, 


de 


—  =  -Y^wt(n)wx(m) 


da[  (k) 


n,m 


p(n,m)  —  Qi(&)  exp  ji  aj(A;)[n  -  /?(m)u i]  } 


The  temporal  and  spatial  windows  are  defined  as 

N-l 

Wt(a)  =  ^>2  wt(n)  exp (— i  an)  for  all  o 

n=0 

M-l 

W*(7)  =  ^2  wx(m)  exp(i  /3  (771)7)  for  all  7 

771=0 


(3) 


for  k  =  ka:kb. 

(4) 

Wt(  0)  =  1, 

(5) 

Wx(0)  =  l, 

(6) 
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and  the  two-dimensional  data  spectrum  as 


N- 1  M- 1 


P(a,7)  =  EE  wt(n)wx(m)  p(n,  m)  exp  [-i  an  -I-  i0(m)j]  for  all  a,y. 


(7) 


n=0  m— 0 


Then,  from  equation  (4),  the  conditionally  optimal  complex  amplitudes  {a:  (k)}kkb ,  for  hypothesized  arrival 
angle  ui ,  satisfy  the  simultaneous  linear  equations 


kb 

Y  a^k)  Wt(ai(k)  -ai(k))  Wa([ai(4)  -aiWjui) 


k=ka 


—  P (an  (k) ,  ai (k)ui ) ,  for  k  =  ka  :kb- 


(8) 


Henceforth,  only  the  case  of  flat  temporal  weighting  is  considered: 

wt(n)  =  for  n  =  0:iV-l. 

Then,  from  equation  (5), 

Wt(a)  =  exp  (  sm(Na/2) 

tK  )  Pl  2  ’N sin(a/2)’ 

from  which 

Wt( 2rt/W)  =  {J  tefcjo}'  W<JV- 
Also,  the  fitting  frequencies  are  taken  as 


(9) 


(10) 


(11) 


fi(k)  =  for  k  =  ka:kb,  ai{k)  =  2nk/N, 

which  leads  to 

Wt(ai(AO-c*i(fc))==Wt  (j^r(k-k)J  =  S(k  -  k). 


(12) 


(13) 


Now,  an  explicit  solution  to  equation  (8)  can  be  written  for  the  conditionally  optimal  amplitudes;  namely, 
Q-tik)  =  P(ai(k),  ai(k)ui)  for  k  =  ka : kb .  (14) 


The  conditionally  minimal  error  then  becomes 


e=^2  wt{n)  wx(m ) 


p(n,m)  -  '^2a1(k)  exp  a\  (k)  [n  -  /3(m)  iti]  j 

TlyTTL  .  k 

=  J2wt(n)  wx(m)\p(n,m)  | 2 

71,771 

-  (k)  Y,Mn)  wx(m)  ]?(n,m)  expjiai(fc)  [n  -  /3(m)  ui]  j . 


p\n,  m) 


71,771 


Pm(ai(k),ai(k)ui)  from  equation  (7) 


(15) 
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Now  the  conditionally  optimal  amplitudes  (equation  (14))  are  substituted  to  obtain 


e  =  ^tot(n)  wx (m)  |p(n,m)|2  -  K  (/c)  |2 

n,m  k 

=  Y^wt(n)  wx(m)  [p(n,m)|2  -  ^|p(ai(fc),ai(fc)ui)| 


=  r(ui) 


(16) 

(17) 


Next,  the  error  e  is  further  minimized  by  choosing  the  planewave  arrival  angle  ui  that  maximizes 
r(u i).  Using  the  definition  of  the  two-dimensional  data  spectrum  P  (equation  (7))  and  the  flat  temporal 
weighting  (equation  (9)),  results  in 


r(ui)  =  ^2 
k 

kb 

=  £ 


IP  ^2wx(m)  P(n,m)  exp  -i^-kn  +  i/3(m)^-ku i 


k=ka 


n,m 

M-l 


wx(m)  q(k,m)ex p 


where 


m= 0 


N- 1 


N 


'.2nk 


N 


(18) 


q(k,m)  =  —  ^  p(n,m)exp(—i2nkn/N )  for  k  =  0:N-l,m  =  0:M-1 


(19) 


n=0 


is  the  temporal  discrete  Fourier  transform  of  the  m-th  element  data.  The  maximization  of  r(m)  by  choice 
of  Hi  is  depicted  in  the  following  sketch: 


r(ui) 


Ui 


Ul 


Best  estimate  of  arrival  angle  is  Ui . 


The  processing  given  by  equation  (18)  has  a  very  plausible  form.  It  says  to  (1)  first  transform  the  time- 
space  data  {p(n,m)}  into  the  frequency-space  domain  {q(k,m)},  and  then  (2)  for  hypothesized  arrival 
angle  U\,  scale  and  phase  shift  the  m-th  element  component  in  frequency  bin  k  by  {2~k/N)j3(m)ui.  But 
this  phase  shift  exactly  compensates  for  that  of  a  single-frequency  planewave  arriving  at  angle  u\\ 


phase  shift (k,m)  =  -2n 


k  x(m)  2irk  n,  . 

«A— — V«m)ui' 


(20) 


Thus,  the  inner  complex  sum  over  m  in  equation  (18)  is  a  coherent  one  for  any  planewave  arriving  at 
angle  tti-  Finally,  the  outer  sum  over  k  in  equation  (18)  is  an  incoherent  sum  over  frequencies  in  the  band 
of  interest.  This  incoherent  sum  is  necessary  because  no  interrelations  have  been  assumed  for  individual 
frequency  components  in  the  planewave  arrivals.  After  the  sum  over  k  is  complete,  then  r(ui)  is  plotted 
for  all  ui  in  the  angular  sector  of  interest,  and  its  maximum  is  located  at  ui. 
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The  complex  sum  on  element  number  m  in  equation  (18)  cannot  be  accomplished  by  a  fast  Fourier 
transform  (FFT)  because,  in  general,  /J(m)  =  x(m)/(cA)  is  not  linear  in  m  for  an  unequally  spaced  line 
array  (but  see  section  2.1.2  for  a  sparse  equispaced  line  array).  This  complex  sum  must  be  carried  out  by 
brute  force.  However,  one  shortcut  available  uses  the  recursion 


exp 


'.27Tfc 

%-jfpym)u  1 


=  exp 


.2n(k  —  1) 
1  N~ 


P(m)u  i 


exp 


.  2tt  .  , 


(21) 


for  each  m  and  u\  to  generate  the  k- values  needed  for  the  exponentials. 


Finding  the  best  spatial  weights  { wx(m )}  is  not  trivial;  they  should  not  simply  be  taken  as  flat  but 
should  reflect  known  element  locations  (x(m)}.  For  example,  wx(m )  might  be  taken  to  be  proportional 
to  the  shaded  area  between  adjacent  element  midpoints,  as  shown  in  the  following  sketch: 


i _ Flat  weighting 


wJm) 


Desired  weighting, 
e.g.,  Hann  or  flat 


x(m  —  1)  x(m)x(m+l) 


After  the  conditionally  best  coefficients  {^(fc)}  are  found  for  a  specified  u\ ,  the  minimal  conditional 
time-space  residual  is,  from  equations  (1)  and  (2), 


p(n,m )  —  p(nA.,x(m))  —  pi(nA,x(m)) 


=  p(n,  m)  -  Oj  (k)  exp 


The  corresponding  FFT  of  this  residual  is,  from  equation  (19), 

1  ^ 

q(k,m)  =  —  ^p(n,m)  exp(-i2Trkn/N) 


n=0 


=  q(k,  m)  -  —  ]T  fij  (k)  exp 


k=ka 


.  27rfc 
-t—Pim) 


n  N- 1 


exp  — z27r(fc  —  k)  n/N ] 


J  n=0 


N6(k  —  k) 


=  q{k,  m)  -  aj(fc)  exp 


.27 rk 

-j— /3(m)ui 


for  k  =  ka:kb. 


The  conditionally  optimal  amplitudes  are,  from  equations  (14),  (9),  (12),  (7),  and  (19), 
oj  (k)  =  P(a1(k),ai(k)ui) 


M- 1 


=  ^  wx(m)  q(k,m)ex p 


771—0 


.271 -k 


for  k  =  ka:kb. 


(22) 


(23) 


(24) 


These  complex  amplitudes  should  be  evaluated  only  after  the  best  arrival  angle  u\  (namely  u\)  has  been 
determined.  Then,  equation  (23)  should  be  evaluated,  using  ui  in  place  of  iq.  That  is,  the  unconditionally 
optimum  amplitudes  are 


ai(fc)  =  P(a1{k),oti(k)iii) 


M- 1 


=  ^  wx{m)  q(k,m)exp 


771=0 


,2nk  „ 


for  k  =  ka: kb, 


(25) 
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and  the  minimal  residual  in  the  frequency-space  domain  is 
q(k,m)  = 


J  q(k,m)  -  fii(fc)exp 

-i^/3(m)ui 

[q(k,m) 

for  k  =  ka:kt,, 
for  k  g  ka:kb . 


Letting  the  spatially  weighted  FFT  at  frequency  bin  k  be  defined  as 
qw(k,m)  =  q(k,m )  wx (m), 

and  its  spatial  autocorrelation  for  frequency  bin  k  be  defined  as 

m 


then,  from  equation  (18),  if  x(m)  =  dm  (equally  spaced  line  array), 
|M— 1  r  t'2 

rM  =  2 


52  ^(fe-m)exp 

fc=fca  |m=0  1  I 

=  E  E  *(M)«p(*^^i«»)- 

*=*„  i=i-M  x  ' 


A  FIM-like2  6  generalization  would  be  to  weight  the  spatial  sum  over  j 
r,(u i)  =51  E  w»W)  ^(*»j)exp  (i-friw i)  » 

j=l-M  '  ' 

where  the  {u>«  0')}^"^  are  separation  weights  (real  and  even  about  j  =  0). 


2.1.2  Sparse  Equispaced  Line  Array 


From  equation  (18)  with  z(m)  =  wx(m)q(k,m)  (holding  A  fixed),  let 
/(«0  =  52  *("») exp[ti/0(m)],  0(m)  =  ^  =  ^«i, 

m=0 

and,  as  an  example  of  a  sparse  equispaced  line  array,  consider 

(x(m)}  =  [x(0)  x(l)  •••  x(M  -  1)] 

=  [0  d  3d  4 d  ^  7d  8 d  9 d  •••]; 

1  missing  2  missing 

then 

{/3(ro)}  =  0  [0  1  3  4  7  8  9  •••],  0  = 

and 


/(i/)  =  z(0)  +  z(l)  exp(iv0)  +z(2)exp(ii'/?3)  +  z(3)  exp(ti/^4) 

+  z(4)  exp(ii/|S7)  +  z(5)  exp(iv/?8)  4-  z(6)  exp(ivj39)  +  ■ 


(26) 

(27) 

(28) 

(29) 

to  obtain  the  modification 

(30) 

(31) 

(32) 

(33) 

(34) 
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Letting 


y-  [z(0)  z(  1)  0  z( 2)  z( 3)  0  0  z( 4)  z(5)  z(6)  •••], 

then 


(35) 


/(i/)  =  y(0)  + 1/(1)  exp(ii//3)  +  2/(3)  exp (ii^3/3)  +  2/(4)  exp(i^4/3) 

+  2/(7) exp(ii/7/3)  +  y(S)  ex.p(iv8(3)  +  2/(9) exp(i^9/3)  +  ••• 

j= o 

=  X^c?')exp 

j=0 


,n  k  d 

,2’nZ a“,j 


(36) 


By  taking  =  jtfr,  this  process  can  be  performed  by  an  iV'-point  FFT.  Sequence  {y(j)}  has  zeros 

in  it,  dictated  by  the  missing  element  locations;  as  a  result,  the  {z(m)}  are  merely  “spaced  out”  by  the 
missing  element  locations. 


2.2  IMPLEMENTATION 


The  theory  outlined  in  the  last  section  was  implemented  in  a  MATLAB  program  and  tested  on  simulated 
data  for  a  variety  of  scenarios. 


2. 2. 1  Coding 


The  first  step  was  to  generate  the  data  that  an  unequally  spaced  line  array  of  sensors  would  produce 
in  the  presence  of  a  planewave  arriving  from  angle  u\  =  sin#i.  First,  the  source  direction,  frequencies, 
and  amplitudes  were  defined: 


ul_source  =  .5;  '/,  sin(source  angle) 

fl_source  =  45:60;  '/,  frequencies  (Hz) 

al=  [0  12345432101234  5];  7.  amplitudes 

The  fitting  band  and  sensor  parameters  were  also  defined: 

M  =  16;  '/.  Number  of  sensors 

wx_m  =  ones(l,M)/M;  '/.  spatial  error  weighting 
ka  =  45; 
kb  =  60; 
i  =  sqrt(-l) ; 

N  =  1024;  '/.  Number  of  time  points 

delta  =  1/1024;  '/,  sampling  increment  (s) 

c  =  1500;  '/,  sound  speed  (m/s) 

tvect  =  delta*(0:N-l)  ;  '/,  vector  of  sample  times 
rand(’ state ’  ,0)  ;  '/,  Same  random  array  each  time 

xarray  =  sort (50*rand(M,  1) )  ;  '/,  random  sensor  positions  (m)  (0,50) 
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The  time-space  data  were  then  computed  as  a  sum  of  planewaves  at  the  source  frequencies  f l.source: 


[xmatrix.tmatrix]  =  meshgrid(xarray,tvect) ; 

pnm  =  zeros (N,M) ; 

for  k  =  1 : length (fl_source) 

pnm  =  pnm  +  al(k)*exp(2*pi*fl_source(k)*. . . 

(tmatrix  -  xmatrix*ul_source/c)*i) ; 

end 


The  next  step  was  to  find  the  value  of  ui  that  maximized  r(uj)  given  by  equation  (18).  First,  a  coarse 
grid  search  was  carried  out  over  L  evenly  spaced  points  in  uj  between  —1  and  1: 


L  =  51; 

ul_l  =  linspace(-l, 1,L) ; 

'/,  Calculate  r(ul)  using  equation  (17) 
beta_m  =  xarray/ (c*delta) ; 

r_l  =  -calcrl(ul_l,  beta_m,  ka,  kb,  N,  wx_m,  pnm); 


A  subroutine  calcrl  was  used  to  compute  the  values  of  r_l  for  an  input  vector  ul_l,  given  the  parameter 
list  beta_m,  ka,  kb,  N,  wx_m,  and  pnm  (see  section  2.2.2).  The  negative  sign  was  used  because  the  function 
calcrl  was  employed  as  the  argument  of  a  function-minimizing  routine  (see  below),  whereas  it  was 
required  to  find  the  function’s  maximum.  The  maximum  value  of  r_l  resulting  from  the  coarse  search 
and  the  corresponding  value  of  ul_l  were  found,  as  well  as  the  ul_l  values  immediately  adjacent: 


[dum.ind]  =  max(r_l); 
if  ind==l 

11  =  1; 

12  =  3; 

elseif  ind==L 

11  =  L  -  2 

12  =  L; 
else 

11  =  ind  -  1; 

12  =  ind  +  1; 

end 


(Maxima  at  the  domain  limits  were  handled  as  special  cases.)  An  accurate  estimate  of  the  peak  was 
obtained  by  using  an  iterative  parabolic  fitting  routine  (fminbnd): 


options  =  optimsetC’TolX’ ,le-10) ; 

best.ul  =  fminbndO  calcrl’  ,ul_l(il)  ,ul_l(i2)  .options  ,beta_m, .. . 
ka,kb,N,wx_m,pnm) ; 

The  fitting  routine  used  the  absolute  tolerance  set  by  the  function  OPTIMSET,  here  equal  to  10-10. 

After  the  best  estimate  of  arrival  angle  (best.ul)  was  found,  the  optimal  amplitudes  {aj (ft)}  were 
then  calculated  from  equation  (14): 
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a_l_bar  =  zeros (1, length (k)) ; 
[beta_m_mx,n_nix]=meshgrid(beta_m,0:  (N  -  1)); 
ik  =  0; 

for  this_k  =  k 

ik  =  ik  +  1;  '/,  index  for  a_l_bar; 
alpha_l_k  =  2*pi*this_k/N; 

Enm  =  exp (( -alpha. l_k*n_mx  +  beta_m_mx*alpha_l_k*best_ul)*i) ; 
Gnm  =  pmn .  *Enm ; 

Fn  =  wx_m*Gnm sum  over  m 
a_l_bar(ik)  =  sum(Fn)/N;  */,  sum  over  n 
end 


2.2.2  Calculation  of  r(ui) 


The  Fourier  transform  of  the  element  data  was  computed  explicitly  for  each  frequency  component  k. 
A  matrix  product  was  used  to  perform  the  sum  over  n  required  in  the  evaluation  of  the  transform  qkm: 

n  =  (0: (N  -  1)).'; 
r_l  =  zeros (size (ul_l)) ; 
for  k  =  ka:kb 

En  =  exp((-2*pi*k*n/N)*i) ; 

'/.  Fourier  transform  at  frequency  k 
'/.  (matrix  product  does  sum  over  n)  : 
qkm  =  (pnm. ’  )*En; 

E_ml  =  exp((2*pi*k*beta_m*ul_l/N)*i) ; 

'/,  matrix  product  does  sum  over  m: 

r_l  =  r_l  +  abs((wx_m.*qkm. ’)*E_ml)  .  "2; 
end  '/,  loop  does  sum  over  k 
r_l  =  -r_l;  !.  for  use  with  MINimizer 

The  negative  of  the  final  result  was  taken  so  that  the  routine  could  be  used  as  the  argument  of  a  minimizer, 
as  explained  previously. 


2.3  EXAMPLES 


2.3.1  Noise-Free  Simulations 


2.3. 1.1  Scenario  1:  Matched  Frequencies.  The  above  code  gave  rise  to  M  =  16  array  ele¬ 
ments  at  the  following  positions:  xarray  =  [0.92518,8.8133,11.557,20.285,22.235,22.823,24.299, 
30.342,  30.772,  36.910,  38.105,  39.597,  41.070,44.565,  46.091,  47.506].  These  positions  are  plotted 
in  the  following  sketch: 


0 


10 


20  30  40  50 

x(m)  position,  m 
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For  this  example,  the  result  of  the  maximization  of  r  by  choice  of  Ui  is  shown  in  figure  1.  In  this  noiseless 
case,  the  resulting  estimate  of  ui  was  0.5  +  5.8  x  10-11,  yielding  the  optimal  amplitudes  {fii(A:)},  shown 
in  figure  2. 

The  optimal  amplitude  errors,  shown  in  table  1,  are  defined  as  ai(k)  —  aj (k).  The  real  parts  of  the 
fitted  amplitude  errors  are  on  the  order  of  the  machine  precision,  and  the  imaginary  parts  are  of  a  size 
comparable  to  the  tolerance  on  the  estimation  of  the  arrival  angle  (10~10),  except  for  the  cases  of  zero 
amplitude,  where  they  are  of  a  size  comparable  to  the  machine  precision.  As  shown  in  the  lower  plot  of 
figure  2,  the  amplitudes  of  the  imaginary  parts  are  proportional  to  the  real  parts. 

For  the  matched-frequency  case,  the  set  of  source  frequencies  was  equal  to  the  set  of  fitting  frequencies. 
In  the  following  four  scenarios  (subsections  2.3. 1.2  to  2.3. 1.5),  cases  where  the  set  of  source  frequencies  is 
not  equal  to  the  set  of  fitting  frequencies  are  considered. 


Figure  1.  The  Result  of  the  Maximization  of  r  by  Choice  of  u i 


Amplitude, 
Real  Part 


k,  or  Frequency,  Hz 


Figure  Z.  Optimal  Amplitudes  for  Noiseless  Example 
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Table  1.  Optimal  Amplitude  Errors  for  Noiseless  Example 


k,  or 

Frequency  (Hz) 

Input 

Amplitude,  ai(fc) 

Optimal  Amplitude  Error 

Real  Part  Imaginary  Part 

45 

0 

1.7  xlO-17  5.2  x  10-15 

46 

1 

2.7  xlO-14  -3.3  x  10-10 

47 

2 

-3.0  x  10-14  -6.8  x  10-10 

48 

3 

-1.7  xlO-14  -1.0  x  10"9 

49 

4 

2.5  xlO-14  -1.4  x  10-9 

50 

5 

9.8  x  10-15  -1.8  x  10-9 

51 

4 

0  -1.5  x  10~9 

52 

3 

2.1  x  10-14  -1.1  x  10-9 

53 

2 

-1.9  xlO-14  -7.7  x  10-1° 

54 

1 

-2.4  x  10~14  -3.9  x  10-1° 

55 

0 

-7.6  xlO-15  9.1  x  10~15 

56 

1 

-1.7  xlO-14  -4.0  x  10_1° 

57 

2 

7.5  x  10-15  -8.2  x  10-1° 

58 

3 

-1.6  xlO-14  -1.2  xlO-9 

59 

4 

-2.9  x  10"14  -1.7  x  10"9 

60 

5 

-3.9  x  10-14  -2.2  x  10“9 

2.3. 1.2  Single  Source  Component  Midway  Between  Fitting  Bins.  In  practice,  the  exact 
source  frequencies  are  not  known  prior  to  conducting  the  single  planewave  fit,  or  the  frequency  bins 
available  in  an  operational  system  may  not  match  those  of  the  source.  To  examine  such  a  situation, 
the  fitting  procedure  was  run  for  a  source  frequency  midway  between  two  of  the  fitting  frequencies,  as 
shown  in  figure  3.  The  input  angle  for  this  case  was  in  =  0.5  (dashed  line  in  the  graph),  and  the  input 
frequency  was  52.5  Hz.  The  fitting  frequencies  were  fk  =  k  =  45, 46, . . .  60.  Although  the  peak  of  the 
function  r(in )  was  found  within  the  requested  tolerance  (10-10),  the  result  differed  from  the  input  angle 
by  3.0  x  10~4.  The  fitted  amplitudes  in  figure  4  show  the  leakage  phenomenon  typically  encountered  with 
Fourier  transform  processing. 


r/maxO)  -  1  q 
(x10-6)_2 

-4 

-6 

-8 


0.499 


“1=sin(91) 
0.4995  0.5 


0.5005 


Figure  3.  Source  Frequency  Midway  Between  Two  Fitting  Frequencies 


2. 3. 1.3  Source  Components  Outside  the  Fitting  Band  and  Misaligned  with  Fitting  Bins. 
To  investigate  the  effect  of  planewave  arrivals  at  frequencies  outside  the  fitting  band,  a  set  of  planewaves 
at  frequencies  above  and  below  the  fitted  frequencies  was  included.  The  resulting  angle  estimate  is  shown 
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in  figure  5.  The  input  angle  for  this  case  was  again  u\  =  0.5  (dashed  line).  The  fitting  frequencies  were 
/*  =  &  =  45, 46, . . .  60.  Although  the  peak  of  the  function  r(ui)  was  found  within  the  requested  tolerance 
(lO-10),  the  result  differed  from  the  input  angle  by  9.5  x  10-4. 

The  fitted  amplitudes  are  shown  in  figure  6.  The  fitted  amplitudes  are  generally  larger  than  the  input 
amplitudes. 


2.3. 1.4  Source  Components  Outside  the  Fitting  Band  and  Aligned  with  the  Fitting  Bins. 
Figure  7  shows  the  effect  of  frequencies  outside  the  fitting  band,  but  this  time  the  actual  and  fitted 
frequencies  are  aligned  over  the  same  equispaced  Fourier  bin  frequencies.  The  angle  is  able  to  be  accurately 
estimated  by  the  fitting  procedure,  as  are  the  amplitudes. 


2. 3. 1.5  Source  Components  Outside  the  Fitting  Band  Are  Misaligned,  While  Source  Com¬ 
ponents  Inside  the  Fitting  Band  Are  Aligned  with  Fitting  Bins.  If  the  frequencies  outside  the 
fitting  band  are  not  aligned  over  the  Fourier  bin  frequencies,  as  in  figure  8,  then  the  optimal  angle  and 
amplitudes  are  not  able  to  be  accurately  estimated  by  the  one-planewave  fitting  procedure. 


Amplitude, 

Real  Part 

0. 


0. 

Amplitude, 
Imaginary  Part 

-0. 


Figure  4 ■  Fitted  Amplitudes  from  One-Planewave  Fitting  Procedure 


u=  sin(0j) 

r/max(r)  -  1  q 
(x10"6)_2 

-4 
-6 
-8 

Figure  5.  Effect  of  Planewave  Arrivals  Outside  the  Fitting  Band 
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Figure  6.  Fitted  Amplitudes  with  Source  Frequencies  Outside  the  Fitting  Band  and 

Misaligned  with  the  Fitting  Bins 


r/max(r)  -  1  *jj. 99 
(x10~6)_2. 

-4- 
-6 
-8i 


w^sinOj) 

0.4995  0.5  0.5005 


40  45  50  55  60  65 


Figure  7.  Results  of  Fitting  with  Source  Components  Outside  the  Fitting  Band  and 

Aligned  with  the  Fitting  Bins 
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(x10~6)_2 

-4 


0.499 


KjSSinCep 

0.4995  0.5  0.5005 


Amplitude, 
Imaginary  Part 


-2\ 


x  Actual 
O  Fitted 


©99 


Q 


- - , - - - _ 

40  45  50  55  60 

k,  or  Frequency,  Hz 


65 


Figure  8.  Results  of  Fitting  with  Source  Components  Outside  the  Fitting  Band  Misaligned, 
While  Source  Components  Inside  the  Fitting  Band  Are  Aligned  with  Fitting  Bins 


It  is  concluded  from  this  set  of  numerical  trials  that  the  fitting  procedure  is  able  to  accurately  estimate 
the  angle  and  amplitudes  of  a  single  multifrequency  planewave,  provided  that  the  Fourier  bin  frequencies 
match  the  frequencies  emitted  by  the  source. 


2.3.2  Comparison  with  Conventional  Beamformer 


In  this  section,  the  one-plaaewave  fitting  algorithm  is  compared  to  conventional  beamforming.  Consider 
the  quantity  r(ui)  (equation  (18)),  which  is  maximized  by  the  choice  of  arrival  angle  u\,  and  let  the 
summand  in  equation  (18)  be  denoted  by  r(u\,k);  that  is, 


*6 

r(ui)=  73 

tefc. 

ki 


M- 1 


73  wx (rn)  q(k,m)  exp 

m= 0 

53  r(wi.fc)> 


.27rk  . 


k=ka 


(37) 

(38) 
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where 


r(ui,k) 


M- 1 

wx(rn)  q(k,m)ex p 

m=0 


.2ttA: 

r— 5(m)ui 


(39) 


The  conventional  beamformer  produces  an  output  at  frequency  index  fc  and  arrival  angle  itj  =  sin(#i ) 
given  by 


1 


Pc(ui,k)  =  j^v'(k)R(k)v(k), 

where  v(fc)  is  the  M-element  steering  vector  at  frequency  index  k  and  at  arrival  angle  6i , 


v(k)  =  exp 
=  exp 


i  2tt — ^  ■  sin(^i)  :r(m) 
c 


'.Stt/c 


N 


(40) 

(41) 

(42) 


and  the  M  x  M  sample  covariance  matrix  R(fc)  =  X(fc)X'(fc),  where  X(fc)  is  the  M-element  vector  of 
the  Fourier  coefficients  at  frequency  fc.  If  equation  (40)  is  compared  with  equation  (18).  it  is  seen  that, 
to  within  a  scaling  factor  N,  the  conventional  beamformer  is  equal  to  the  summand  of  the  function  r(u\) 
with  flat  spatial  weighting. 


In  the  following  examples,  the  “matched  frequencies”  case  of  the  previous  section  is  used;  that  is,  the 
planewave  is  composed  of  frequencies  corresponding  to  k  =  fk  =  45:60,  and  the  arrival  angle  corresponds 
to  ui  =  0.5. 


2.3.Z.1  Noise-Free  Case.  Figure  9  shows  the  function  r(ui,k)  (equation  (39))  on  the  left  and 
the  conventionally  beamformed  output  (equation  (40))  on  the  right.  The  darker  regions  denote  high 
amplitude,  and  the  lighter  regions  denote  low'  amplitude.  The  set  of  input  frequencies  is  shown  as  a  set 
of  white  dots  at  m  =  0.5.  The  intensity  image  was  produced  using  101ogr(u1,k).  For  the  conventional 
beamformer  showm  at  the  right,  temporal  FFTs  with  a  length  equal  to  the  number  of  time  samples  for 
each  sensor  (1024  points)  were  used  to  estimate  the  covariance  matrix.  Both  plots  have  the  same  k,  ui, 
and  gray  intensity  scales.  As  can  be  seen,  the  two  plots  sire  virtually  identical. 
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Figure  9.  Comparison  of  Function  r(ui,k)  ( Equation  (39))  (Left)  and  Conventionally 

Beamformed  Output  (Right) 
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Figure  10  shows  r(ui)  and  the  sum  over  k  of  the  conventional  beamformer  output,  both  plotted  as 
functions  of  u\  on  a  linear  ordinate  scale.  The  function  r(u i )  corresponds  to  the  conventional  beamformer 
output  summed  over  the  fitting  frequencies.  The  resulting  peak  is  at  the  correct  value  of  ui;  namely, 
uj  =  0.5.  The  two  plots  are  identical  and  hence  superimposed. 

Figure  11  shows  the  strengths  of  the  fitted  frequency  components  and  the  amplitudes  of  the  conven¬ 
tional  beamformer  output  plotted  against  the  frequency  bin  number  at  u\  =  0.5.  The  top  plot  shows  the 
actual  strengths,  and  the  bottom  plot  shows  (1)  the  differences  between  the  one-planewave  fit  and  the 
input  strengths  (circles),  and  (2)  the  differences  between  the  conventional  beamformer  output  and  the 
input  strengths  (triangles).  There  is  a  small  error  between  the  conventional  beamformer  output  and  the 
input  strengths,  which  is  not  evident  in  the  one-planewave  fit. 
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£ 
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Figure  10.  The  Quantity  r(u\)  and  the  Sum  over  k  of  the  Conventional  Beamformer 

Output,  Both  Plotted  as  Functions  uf  u\ 
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Figure  11.  Strengths  of  the  Fitted  Frequency  Components  and  the  Amplitudes  of  the 
Conventional  Beamformer  Output  Plotted  Against  Frequency  Bin  at  Ui  =  0.5 
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2.3. 2. 2  Case  with  Noise.  Complex  normally-distributed  random  noise  is  added  to  the  simulated 
pressure  data  in  the  following  manner: 


snr  =  -10; 

Noise  =  randn(size(pnm))  +  i*randn(size(pnm)) ; 
As  =  mean (std (pnm) ) ; 

An  =  mean (std (Noise)) ; 

NoiseScaleFactor  =  As  /  (An*10~(snr/20)) ; 
pnm  =  pnm  +  Noise  *  NoiseScaleFactor; 


The  quantity  snr  defines  the  desired  signal-to-noise  ratio,  which  is  equal  to  —10  dB  in  this  case.  The 
matrix  Noise  is  the  same  size  as  the  data  matrix  and  has  independent  real  and  imaginary  parts.  The 
quantity  std  (pnm)  is  an  M-element  vector  of  standard  deviations — one  for  each  array  element.  The 
quantity  As  is  the  mean  of  these  standard  deviations:  a  similar  measure  is  made  for  the  noise  amplitude 
An.  The  noise  matrix  is  multiplied  by  the  scalar  NoiseScaleFactor  to  obtain  the  desired  SNR,  and  this 
quantity  is  added  to  the  data  sample  matrix. 

Figure  12  shows  the  function  r(ui,k)  (equation  (39))  on  the  left  and  the  conventionally  beamformed 
output  (equation  (40))  on  the  right  for  this  case,  with  an  SNR  of  -10  dB.  For  the  conventional  beamformer 
shown  on  the  right,  temporal  FFTs  with  a  length  equal  to  the  number  of  time  samples  for  each  sensor  (1024 
points)  were  used  to  estimate  the  covariance  matrix.  The  two  plots  are  very  similar.  Further  similarities 
between  the  one-planewave  fitting  algorithm  and  conventional  beamforming  are  identified  below. 


Figure  12.  The  Function  r(u\,k)  (Equation  (39))  (left)  and  the  Conventionally 
Beamformed  Output  (Equation  (40))  (right)  for  SNR  of  —10  dB 


Figure  13  displays  r(u i)  and  the  sum  over  k  of  the  conventional  beamformer  output,  both  plotted  as 
functions  of  ui  on  a  linear  ordinate  scale.  The  two  lines  are,  as  previously  shown,  virtually  identical. 

Figure  14  shows  the  frequency  strengths  plotted  against  the  frequency  bin  number  at  m  =  0.5.  As 
seen  earlier  in  the  noise-free  case,  there  is  no  significant  difference  between  the  one-planewave  fit  (dots) 
and  the  conventional  beamformer  (triangles). 

Figure  15  compares  the  outputs  of  the  conventional  beamformer  and  the  one-planewave  fit  for  the 
range  of  input  SNRs  shown  at  the  left  of  each  row  of  plots.  Ten  realizations  are  shown  for  each  SNR  value. 
The  axis  scale  for  each  column  is  shown  at  the  bottom  and  is  constant  within  each  column,  except  for  the 
second  one.  The  first  column  show’s  the  conventional  beamformer  output  plotted  as  a  surface  in  decibels; 
the  gray  scale  values  are  indicated  to  the  left  of  the  top  plot  and  are  constant  throughout  column  1.  The 
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input  frequencies  are  shown  as  the  white  dots,  and  the  peak  value  of  the  beamformer  output  is  shown  by 
the  white  cross  in  each  plot.  Column  2  shows  r(«i)  and  the  sum  over  k  of  the  conventional  beamformer 
output,  both  plotted  as  functions  of  ui  on  a  linear  ordinate  scale.  Solid  lines  in  this  column  are  the 
one-planewave  fits;  broken  lines  show  the  sum  over  k  of  the  conventional  beamformer  outputs,  but  these 
axe  superimposed  on  the  solid  lines  and  are  hence  invisible.  Column  3  shows  the  frequency  strengths 
for  the  one-planewave  fit  (circles)  and  the  conventional  beamformer  (triangles)  at  Ui  =  0.5.  The  set  of 
amplitudes  listed  in  table  1  on  page  12  was  used.  There  is  no  significant  difference  between  the  outputs 
of  the  one-planewave  fit  and  the  conventional  beamformer. 


ui  =  sin(0i) 


Figure  13.  The  Quantity  r(uj)  and  the  Sum.  over  k  of  the  Conventional  Beamformer 

Output,  Both  Plotted  as  Functions  of  ui 


k,  or  frequency  in  Hz 


Figure  14 •  Frequency  Strengths  Plotted  Against  Frequency  Bin  Number  at  u\  =  0.5  for  the 

Case  with  Noise 
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3.  TWO-PLANEWAVE  FIT  FOR 
UNEQUALLY-SPACED  LINE  ARRAYS 


3.1  THEORY 


The  discrete  data  available  are  again  {p(n,m)}  for  n  =  0:N—  1  and  m  =  0:M-1,  and  the  fitting 
frequencies  are  taken  as 


k  k 

h  T  NA 


for  k  =  ka:kb. 


(43) 


The  range  (ka,  kb)  represents  the  common  frequencies  for  the  two  planewave  fits.  (It  is  possible  to  generalize 
to  the  case  where  only  some  of  the  fitting  frequencies  {/*}  are  common  to  the  two  planewaves;  the 
corresponding  complex  amplitudes  {ai(A:)}  and  {02 (A:)}  would  then  be  set  to  zero  in  the  appropriate 
nonoverlapping  frequency  bins.)  Also  defined  are 


_  2 nk  _  x(m) 

a*  =  -jj-,  «m)  =  -3- 


for  m  =  0 :  M  - 1. 


(44) 


If  two  planewaves  arrive  from  hypothesized  angles  u\  and  u2,  the  observed  complex  pressure  field  over 
observation  interval  T  =  NA  is  modeled  as 


kb  kff 

P2 (t,x)  =  53  oi(*)exp[i2fff  (f-  fui)]  +  53  fl2(A:)exp[i27r^  (t-  §u2)j, 


(45) 


k=ka 


Jfc=Jtn 


where  the  hypothesized  amplitudes  {ai(A:)}  and  (a2(/c)}  are  complex.  The  two  (broadband)  sources  are 
assumed  to  be  uncorrelated  with  each  other  (see  section  4  for  a  correlated  case).  Therefore,  the  sampled 
pressure  field  is  modeled  as 


p2(nA,x(m))  =  53  ai(A:)exp[ia*(n- j8(m)tti)j 
k=ka 
kb 

53  a2(fc)  exp  [ta*  (n  -  0(m)u2)]  ■ 


(46) 


k=ka 


A  total  average  squared  error  between  the  actual  received  data  and  the  model  (equation  (46))  is  defined 


as 


N- 1  M—l 


1  x  *  »  *  1 2 

e  =  —  53  13  w*(m)  P(.nAMm))  -ft(nA,®(m))|  , 


(47) 


n=0  m= 0 


where  a  flat  temporal  weighting  wt(n)  =  1/jV  for  n  =  0:N  —  1  is  used,  and  a  spatial  error- weighting 
function  wx(m),  which  is  real  and  positive,  is  allowed.  Then  the  error  e  can  be  expressed  in  the  form 


n,m 


e  =  —  ^wx{m)  p(n,m)  -  53  (fc) exp [ia* (n  -  0(m)«i)j 

k=ka 

kb  ; 

-  53  a2(Jk)exp^'aifc(n-)3(m)«2)j 


k=ka 


=  \d(n,m)\2  =  Y^v*(rn)jjYl\d(n’m')\2- 


(48) 

(49) 
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Letting 


N-l 


D(l,m )  =  E  exp(— i2-nnl/N)  d(n,m )  for  l  =  0:N  —  1, 


n= 0 


then 


N- 1 


JV-1 


E  EexP 


iV2 


i=0 


i=0  n,n 
N-l 


-i2n(n  —  n) 


d(n,m )  d*(n,m) 


=  ^  E  id(n>m)i2- 


n=0 


which  is  a  discrete  form  of  Parseval’s  Theorem.  The  error  e  follows  from  equation  (49) 
e  =  E^(m)^El^m>l2- 

m  / 

However,  from  equations  (49),  (50)  and  (44), 

jyD(l,  m)  =-^Eexp (-*27rn//iV){p(n,m) 

-E«i(*)«P  io:/b(n  -  /3(m)ux)j  —  Ea2(^)exp  ia:jfc(n  - /?(m)u2) 
k  k 

—q(l,  m)  —  ai(l)  exp[— iaj/?(m)ui]  -  02(1)  exp[— iaj/?(m)it2], 
where  the  temporal  discrete  Fourier  transforms  have  been  defined  as 
1  JV_1 

q(l,m)  =  —  E  exp  (— i2irnl/N)p(n,m), 


n= 0 


for  l  =  0 :  iV — 1,  m  =  0:M  —  1.  Therefore,  equation  (52)  becomes 
e  =  E^*(m)E  |^> m)  -  ai(0  exp  [-*o:i/3(TO)wi] 

m  ( 

2 

-  a2(Z)  exp [-iQ(/3(m)u2]  . 

To  minimize  e,  consider 

=  —  E  -ai(fc)exp[-iafc/3(m)ui] 

-  a2(fc)  exp[— iajtj0(m)u2]  j  exp[iafc/3(m)ui] 
=  -E«-xH  g(A;, m)exp[ia!jb/3(m)ui]  +  Qi(fe) E ^xN 

m  m 

+  a2(fc)  E  wx{m)  exp[iak/3(m)(ui  -  it2)] 

m 

-  -X(k,akui)  +  ai(k)  +  a2(k)Wx(ak(u  1  -  u2))  for  k  =  ka:kb, 


where 


X{k,  7)  =  ^2wx(m)  q(k,m)  exp[i7/?(m)]  for  k  =  ka'.h,  all  7, 

m 

^(7)  =  y^  wx(m)  expfi7/3(m)]  for  all  7, 

m 

and  Wj  (0)  =  ^  wx  (m)  =  1  is  used.  Similarly, 

m 

9e 


(57) 

(58) 


dn2{k) 


=  -X(k,  aku2 )  +  fli  (A)  Ws*  (a*  (ui  -  u2))  +  a2(k) 


(59) 


for  k  =  ka:  kb.  Minimization  of  error  e  results  in  simultaneous  equations  for  the  conditionally  best-fitting 
coefficients  {aj  (A)}  and  (o^A)},  namely, 


fii ( k)  +  W*  a2 (A)  =X(k,akui)\ 
WZ  ax{k)  +  a2(k)  =  X(k,aku2)  j 


>  for  k  =  ka:kb, 


(60) 


and  Wk  =  Wx  (ak(ui  —  u2)).  The  solutions  for  the  conditionally  optimal  complex  amplitudes  are,  for 

k  —  fcfl  *  kb  j 


,1S  X(A,a*iii)  -Wk  X(k,aku2) 
~l[  )  ~  1  -  |W*|a 

,IN  X(k,aku2) -WZ  X(k,aku  1) 

-2()"  i-F*l2 


(61) 


Decoupling  in  frequency  index  k  is  due  to  the  flat  temporal  weighting  {u>t(n)}  and  to  the  fact  that 
frequency  increment  A /  =  1/T.  The  conditionally  minimal  value  of  e  is  now,  from  equation  (55), 


e  =  X) Wx  (m)  { 9(*’ m)  -  £1  ( k )  exP  [-iakP{m)ui] 

m  k 

-  02(A)  exp  [-iak/3(m)u2]  |g*(A,  to) 

=  5Z^i(m)  |g(A, m)|2  -  ^a1(A)X*(A,afc'Ui)  -  a2(A)X*(A,a*u2) 


m,fc  fc 

=  ^2wx(m)  |g(A,m)|2  -r(m,u2), 

rriyk 


(62) 


where  the  solutions  for  a^A)  and  a2(A)  from  equation  (61)  have  been  used,  as  well  as  the  fact  that  the 
initial  double  sum  over  m  and  A  is  independent  of  in  and  u2,  and 


r(ui,u2)  =  ^ 


|X(A,^u1)|2  +  |X(A,^u2)j2 
^  -2»{w,(2j^(«1  -u2))X*(A,^Ul)X(A,2^U2)} 


k=ka 


Wx(2#(u  1  -ua)) 


(63) 
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where  9?  denotes  the  real  part,  and,  from  equations  57,  58,  and  54, 


M- 1 


X(k,  - fru )  =  X)  wx(m)  q(k,m)ex p 

m= 0 
Af-1 

=  5Z  wx(m)exp 


\.2nk  x(m)  ' 


m— 0 

iV-l 


.2ttA;  x(m) 

j - 1 

iV  cA 

,  (W*(0)  =  1), 


N  cA 

1  1'~± 

q(k,  m)  =  —  ^2  P(n > m )  exp(—i2nkn/N)  for  m  =  0 : M- 1 


n=0 


(64) 

(65) 

(66) 


for  A:  =  0:1V-1.  However,  only  the  values  for  k  =  ka-h  are  needed  in  equation  (63).  Equation  (66) 
consists  of  M  JV-point  temporal  discrete  Fourier  transforms. 


To  further  minimize  the  conditional  error  e  in  equation  (62),  the  quantity  r(ui,  u2)  given  by  equa¬ 
tion  (63)  must  be  maximized  by  choice  of  both  arrival  angles  ui  and  1x2  >  which  may  be  taken  as  ui  <  u2 
without  loss  of  generality.  The  location  of  the  maximum  of  r(ui,u2)  is  denoted  by  ui,U2- 

The  factor  in  equation  (64)  is  exactly  the  phase  compensation  required  at  frequency 

and  array  element  m,  to  coherently  “line  up”  all  components  arriving  at  angle  u  (u  =  0  is  broadside). 
Thus,  equation  (64)  is  a  “coherent”  sum  carried  out  prior  to  the  “incoherent”  sum  over  k  (frequency)  in 
equation  (63).  The  sum  for  spatial  window  Wx  in  equation  (65)  cannot  be  carried  out  a  priori  in  closed 
form  because  of  the  irregular  element  locations  (x(m)};  also,  ^(7)  is  complex  and  will  remain  so  for 
general  (x(m)}  locations. 

From  equation  63,  it  is  noted  that  for  given  ui,  u2,  and  k,  only  three  complex  quantities  need  to  be 
computed.  However,  they  must  be  combined  according  to  equation  (63),  yielding  a  purely  real  quantity, 
which  is  then  summed  over  the  frequency  band  of  interest,  k  =  ka:kt,.  Equations  (64)  and  (65)  cannot  be 
carried  out  with  the  FFT  because  the  element  locations  {x(m)}  are  unevenly  spaced. 

A  measure  of  the  total  power  in  arrival  1  is  available,  by  reference  to  equation  (45) ,  as 
P1=  £  |a!(fc)|2,  (67) 

k—ka 

where  {di  (fc)}  are  the  optimal  coefficients  obtained  from  equation  (61),  using  angle  values  Ui  and  112,  after 
the  best  angles  (namely,  Ui  and  u2)  have  been  determined  from  the  maximization  of  r(u1;u2)  given  by 
equation  (63).  If  Pi  is  large,  then  source  1  could  be  coherently  subtracted  from  the  input  data  p(n,m),  as 
indicated  in  equation  (48).  If  Pi  is  small,  then  the  coefficients  {oi(fc)}  could  be  discarded  and  considered 
as  noise. 


With  the  conditionally  optimal  coefficients  {^(k)}  and  (a2(A:)}  for  specified  u\  and  U2,  the  condition¬ 
ally  minimal  time-space  residual  is,  using  equation  (46), 


p(n,  m)  =  p(nA,  x(m))  —  p2{nA,  x(m)) 


=  p{n,m)  -  ^Pa^A)  exp 


p{m)u2) 


(68) 
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The  unconditionally  minimal  time-space  residual  is 


p(n, m)  =  p(n, m)  ~y~] a\(k) exp 

k 

— ^2  0s)  exp 

k 


i~N~  (n-'9(m)*0 
i^(n  -  /3(m)u2) 


(69) 


for  n  =  0 :  TV  —  1,  m  =  —  The  corresponding  unconditionally  minimal  residual  in  the  frequency- 

space  domain  is,  using  equation  (66), 


q(k,m )  =  q(k,m)  —  &i(k)  exp 


.27 rk 
l~N 


P(m)u  i 


277  fc 

—  5,2 (fc)  exp  -i—^-P(m)u2  for  k  =  &a :&&,  m  =  0:M  —  1. 


(70) 


These  residuals  are  devoid  of  the  two  strongest  planewave  arrivals  and  can  now  be  processed  further  to 
detect  additional  (weak)  arrivals.  Without  this  coherent  subtraction,  the  two  strong  arrivals  could  override 
a  weak  arrival  in  close  angular  proximity. 


3.2  IMPLEMENTATION 


The  first  step  of  the  two-planewave  fit  is  to  find  the  pair  of  values,  u\  and  h2,  that  maximize  the 
function  r(ui,u2)  given  by  equation  (63).  Three-dimensional  arrays  are  used  to  store  the  quantities 
X(k,2nkui/N),  X{k,  2irku2/N),  and  Wx(2nk(ui  -  u2)/N),  each  being  considered  as  a  function  of  uj, 
u2,  and  k.  The  array  geometry,  sampling  interval,  planewave  amplitudes  and  directions,  etc.,  are  defined 
as  follows: 


fl .source  =  45:60; 

al  =  [0  12345432101234  5]; 

ka  =  45; 
kb  =  60; 

f2_source  =  45:60; 
a2  =  5  -  al; 

ul.source  =  0;  '/,  sin(source  angle) 

u2_source  =  .5;  7,  sin  (source  angle) 
i  =  sqrt(-l) ; 

N  =  1024;  */.  Number  of  time  points 

delta  =  1/1024;  V,  sampling  increment  (s) 
c  =  1500;  V,  sound  speed  (m/s) 

tvect  =  0:delta:(N  -  l)*delta;  '/.  vector  of  sample  times 
xarray  =  linspace(0,50,ll) . ’ ; 

M  =  length  (xarray) ;  '/,  Number  of  sensors 
wx_m  =  ones(l,M)/M;  7,  spatial  error  weighting 


The  time-space  data  are  generated  by  summing  the  frequency  components  of  the  two  sources: 

[xmatrix ,  tmatrix]  =  meshgr  id  (xarray,  tvect) ; 
pnm  =  zeros (N,M) ; 

argl  =  tmatrix  -  xmatrix*ul_source/c; 
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arg2  =  tmatrix  -  xmatrix*u2_source/c; 

K  =  length(f l_source) ; 
for  k  =  1 : K 

pnm  =  pnm  +  al(k)*exp(2*pi*f  l_source(k)*argl*i)  +  ... 
a2  (k) *exp (2*pi*f 2_source (k) *arg2*i) ; 

end 


The  Fourier  transform  q(k,m )  given  by  equation  (66)  (an  M  x  K  matrix)  and  the  phase  factor  <j>  = 
2-7T  k  x(m) /(NcA)  (also  M  x  K)  used  in  equations  (64)  and  (65)  are  calculated  as  follows: 


*/.  Calculate  Fourier  transform  and  phase  matrix  for  use  in  the 
'/,  calculation  of  r: 
k_vt  =  ka:kb; 
n_vt  =  0 : N— 1 ; 

[k_mx,n_mx]  =  meshgrid(k_vt,n_vt) ; 

E  =  exp((-2*pi*k_mx.*n_mx/N)*i) ; 

7,  Fourier  transform:  sum  over  n  done  with  matrix  product: 
q  =  (1/N)*pnm. J*E; 

V.  Phase  factor  used  in  argument  of  exponential 

y,  in  calculation  of  X  and  W: 

phi  =  2*pi*xarray*k_vt/(N*c*delta) ; 


The  vector  xarray  is  of  size  Mx  1,  and  the  vector  k_vt  is  of  size  lxK.  The  quantity  Wx  (2  r  k  (u\  -u2)/N) 
is  calculated  using  two  nested  f  or-loops  and  stored  in  an  Nu  x  Nu  x  K  array,  while  an  intermediate  array 
Xku  stores  values  of  X(k,  2irkujN)  in  a  K  x  Nu  array,  where  u~- 1  +  <5U[0:(1VU  -  1)],  Su  =  2/(Nu  -  1), 
and  Nu  is  the  number  of  values  used  to  span  the  range  of  u: 


Nuspan  =  40; 

u_span  =  linspace(-l,l, Nuspan); 

Xku  =  zeros (K, Nuspan)  ; 

Wkuu  =  zer  os  (Nuspan,  Nuspan,  K)  ; 
iu2  =  0 ; 
for  this_u2  =  u_span 
iu2  =  iu2  +  1; 

'/,  Matrix  product  does  sum  over  m: 

Xku(:,iu2)  =  (wx_m*(q.*exp(phi*this_u2*i))) .  ’ ; 
iul  =  0; 

for  this_ul  =  u_span 
iul  =  iul  +  1; 

Emkuu  =  exp(phi*(this_ul  -  this_u2)*i); 

V,  Matrix  product  does  sum  over  m: 

Wkuu (iu2, iul, :)  =  wx_m*Emkuu; 
end 
end 


The  K  x  Nu  array  Xku  is  now  copied  into  two  (three-dimensional)  arrays  that  have  the  same  dimensions 
as  Wkuu.  These  dimensions  correspond  to  ui  varying  across  the  columns  and  u2  varying  down  the  rows. 
In  all  cases,  k  varies  along  the  third  dimension.  This  copying  is  carried  out  as  follows: 
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Xkul  =  permute (repmat (Xku, [1  1  Nuspan] ), [3,2, 1] ) ; 
Xku2  =  permute (repmat (Xku ,[11  Nuspan] ) , [2 , 3 , 1] ) ; 


The  matrices  Wkuu,  Xkul,  and  Xku2  are  now  used  to  compute  the  addend  in  equation  (63),  with  the  result 
summed  over  k,  corresponding  to  dimension  number  3: 


top  =  abs (Xkul). “2  +  abs (Xku2) . ~2  ... 

-  2*real (Wkuu. *0011  j  (Xkul)  .*Xku2) ; 
bot  =  1  -  abs (Wkuu) . ‘2; 

rulu2  =  sum  (top. /bot,  3) ;  '/.  sum  over  the  third  index,  which  is  k 


The  matrix  rulu2  is  of  size  Nu  x  Nu.  The  location  of  the  maximum  value  of  this  matrix  gives  an  initial 
estimate  of  u\  and  U2,  as  follows: 


[dum.jmax]  =  max(max(rulu2))  ; 
[dum,  imax]  =  max(rulu2( : ,  jmax) ) ; 
hatul  =  u_span(jmax) ; 
hatu2  =  u_span(imax) ; 


Figure  16  shows  calculations  of  the  quantity  r(u\ ,  112)  for  different  values  of  the  actual  source  angles  us\ 
and  uS2-  Contours  are  plotted  at  the  set  of  values  shown  at  the  top  right  of  the  plots,  with  the  actual  source 
angles  us  3  and  us 2  indicated  by  the  white  dot  in  each  plot.  The  function  r  (7x3, 7x2)  is  symmetric  under 
exchange  of  u\  and  u2,  as  expected  by  equation  (63).  Equation  (63)  is  undefined  at  ux  =  u2;  therefore, 
calculating  r (7x3,112)  for  7x3  =  7x2  has  been  avoided.  (See  the  appendix  for  an  alternative  formulation 
without  this  singularity.)  For  the  cases  where  the  source  angles  were  equal  (plots  at  the  top  and  right 
edge  of  the  figure),  the  r (7x3, tx2)  surfaces  have  no  unique  maxima  and  correspond  to  a  single  planewave 
source. 

Initial  estimates  of  txoi  and  U02  were  obtained  by  finding  the  location  of  the  peak  of  rulu2;  now  those 
estimates  will  be  refined  using  the  iterative  procedure  fminsearch,  and  the  function  calcrulu2,  which 
returns  a  scalar  r(u3,7x2)  for  a  given  input  vector  [txi,  7x2). 


3.3  EXAMPLE  OF  NOISE-FREE  SIMULATIONS 


The  sketch  below  shows  a  noise-free  simulation  that  was  run  for  an  array  of  20  randomly  distributed 
sensors  placed  in  a  line  at  the  following  positions  (meters):  0,  2.7149,  6.3877,  8.2533,  9.5938,  12.3172, 
13.2960,  13.9596, 14.6925,  23.4220,  26.1543,  33.2634,  37.1476,  40.0923,  44.2354,  44.4607,  45.1650,  48.3877’ 
49.5603,  and  50. 
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Two  sources  were  simulated,  one  at  us\  =  —0.01,  and  the  other  at  uS2  =  +0.01,  i.e.,  0.57°  either  side 
of  broadside.  Both  emitted  a  set  of  pianewaves  at  frequencies  between  45  and  60  Hz,  with  a  spacing  of 
1  Hz.  The  amplitudes  of  the  pianewaves  are  as  shown  in  figure  17. 

Figure  IS  shows  contours  of  the  quantity  r(ui,tt2).  The  plot  on  the  left  shows  r  for  all  possible  values 
of  ui  and  u%,  and  the  plot  on  the  right  presents  a  magnified  view  of  a  small  region  around  the  true  angles 
of  the  sources.  The  angles  found  by  the  search  procedure  are  indicated  on  the  magnified  plot,  with  the 
point  corresponding  to  these  angles  laying  on  a  very  shallow  hyperbolically  (banana-)  shaped  ridge  that 
asymptotically  approaches  the  lines  ui  =  0  and  no  =  0.  The  entire  surface  is  symmetric  under  exchange 
of  ui  and  uo,  as  expected. 


Source 

Amplitude 


Figure  17.  Amplitudes  of  the  Two  Pianewaves  Used  for  the  Two-Planewave  Fit 
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Figure  18.  Quantity  r(ui,«2)  for  the  Two-Planewave  Fit 
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The  previous  result  is  compared  with  conventional  beamforming  results  in  figure  19.  The  top  plot 
provides  contours  of  the  same  two-planewave  fit  that  was  shown  in  the  previous  plot,  together  with  the 
zoomed-in  region  to  the  right,  close  to  the  true  source  values.  The  conventionally  beamformed  result  is 
shown  in  the  lower-left  plot,  with  a  zoomed-in  version  to  the  right,  surrounding  the  true  source  angles 
—0.01  and  +0.01.  The  conventional  beamformer  output  was  found  for  a  frequency  of  50  Hz,  which  was 
one  of  the  frequencies  corresponding  to  the  largest  amplitude  emitted  by  the  sources  (see  the  source 
amplitude  plot  in  figure  17).  There  is  no  evidence  at  all  of  the  angles  corresponding  to  the  two  sources  in 
the  conventionally  beamformed  result. 
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Figure  19. 


The  r(ui ,  u-z)  Surface  Compared  to  Conventionally  Beamformed  Results 
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Figure  20  shows  the  same  comparison,  but  this  time,  the  sources  have  been  more  widely  separated  in 
angle,  corresponding  to  us i,u*2  =  ±0.1715.  These  values  were  found  by  adjusting  the  angular  separation 
until  a  double  peak  began  to  be  observed  in  the  conventional  beamformer  output.  The  lower-right  plot 
shows  that  even  though  the  conventional  beamformer  was  able  to  produce  a  double  peak  in  this  instance, 
the  peaks  did  not  occur  at  the  correct  values  of  u,  which  are  indicated  by  the  broken  lines. 
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Figure  20.  The  r(u  1,112)  Surface  Compared  to  Conventionally  Beamformed  Results  for 

Separated  Sources 
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Figure  21  shows  the  results  for  very  widely  separated  sources  and  for  a  single  source.  Case  (a)  illustrates 
the  presence  of  two  sources  (at  ue\,ua2  =  ±0.75)  and  case  (b)  illustrates  the  presence  of  only  one 
source  (at  usl  =  —0.75).  The  top  plots  are  the  two  two-planewave  fits,  and  the  bottom  plots  are  the 
conventionally  beamformed  results.  The  two-planewave  fit  ((a)  top)  was  able  to  accurately  determine  the 
angles  corresponding  to  the  two  plane  waves.  Although  the  conventional  beamformer  ((a)  bottom)  shows 
two  distinct  peaks,  the  mutual  interference  between  the  two  planewaves  gave  rise  to  peaks  that  were 
closer  together  in  angle  than  the  true  sources.  The  same  conventional  beamformer  was  able  to  accurately 
determine  the  direction  of  the  single  planewave  ((b)  bottom). 


(a) 

Two-planewave  fit  and 
conventionally  beamformed  result 
for  widely  separated  sources 


(b) 


Two-planewave  fit  and 
conventionally  beamformed  result 
for  a  single  source 
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Figure  21.  Two-Planewave  Fits  and  Conventional  Beamforming  Outputs  for  (a)  Two 
Widely  Separated  Sources  and  (b)  A  Single  Source 
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4.  TWO-PLANEWAVE  FIT  TO  MULTIPATH  ARRIVALS 
FOR  UNEQUALLY  SPACED  LINE  ARRAYS 

The  discrete  data  available  are  {p(n,m)}  for  n  =  0:-/V-l  and  m  =  0:M-1.  Here,  p(n,m )  = 
p(nA,x(m)),  where  p(t,  x)  is  the  continuous  received  pressure  field  at  time  t  and  location  x,  while  A 
is  the  time-sampling  increment  and  {ar(m)}  are  the  sensor  positions. 

A  planewave  arriving  from  angle  u\  =  sin  $1  is  modeled  according  to  the  following  summation: 


ki 


Pi(t,x)  =  ^2  exP  -  fui)]  . 


(71) 


k=ka 


A  multipath  arrival  with  relative  gain  G  and  delay  r  from  arrival  angle  ii2  is  modeled  as 


kb 


Pi(t, x)  =  G  di(k)  exp  [i2ir±(t  -  r  -  fu2)] , 


(72) 


k=ka 


where  the  same  set  of  complex  amplitudes  (ai(A:)}  is  utilized.  Gain  G  is  real  and  independent  of  frequency 
in  the  band  of  interest.  As  a  result,  the  modeled  samples  at  times  {nA}  and  locations  {x(m)}  comprise 
a  two-planewave  fit: 


P2{nA,x(m))  =  (A:)  exp 


it 


i27r—  nA|  ^  exp 


+Gexp 


k  x (m) 

-,2*t -T' “■ 


..  k  k  x(m)  1 

-t27T-T  - 1^-——  u2 


The  dimensionless  parameters  are  defined  as 


_  2l(k  at  \  -  x(m ) 

“*  =  -JP  ^  =  IF' 


(73) 

(74) 

(75) 


The  modeled  data  samples  can  then  be  written  as 

p2(nA,i(m))  =  ^ai(fc)exp(iafcn)  jexp[— iakP(m)ui]  +  Gexp[-iafcA  -  iak0(rn)u 2]|.  (76) 


An  average  squared  error  between  the  data  and  the  fit  (equation  (76))  is  defined  as 
e  =  ^]>^iUx(m)jp(nA,z(m))  -^(nA^m)) j  , 


(77) 


where  the  temporal  weighting  is  fiat,  and  {wx(m)}  are  real  positive  spatial  error  weights,  which  sum  to  1. 
The  error  can  be  expressed  as 


p(n,  m)  -  ^2  °i (k)  exp(iafcn)  &(fc,m;ui,u2,G,  A) 


where  b(k,m)  =  b(k,m;ui,U2,G,X)  is  given  by 

b(k,m-,ui,U2,G,X)  =  exp[-ia*/?(m)ui]  +  Gexp[-iafcA  -  iafc/?(m)u2] . 


(78) 


(79) 
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Defining  difference 

d(n,  m )  =  p(n,  m)  ~Y  oi  ( k )  exp(ia*n)  b(k,  m)  (80) 

k 

then  yields,  from  equations  (78)  and  (51), 

e  =  Y,Wx^jjH  M(n»m)i2  =  5^I-D(*»m)l2-  (81) 

m  n  ml 


Now,  from  equations  (80)  and  (50), 

1  \  1  ^r— \  (  .2nnl\  . 

-D(l,m)  =  ^  2^  exp  l-i—  1  d(n,m) 

n=D  '  ' 


p(n,m)  -  Y  aiW  exp 


1  /  ,27TTll\ 

=  ivEexP^iv-J 

n  v  ' 

=  q(l,m)  -  ai(l)  b(l,m), 

where  the  {q(f,m)}  are  the  temporal  discrete  Fourier  transforms 


27r&n 

“iv“ 


6(fc,  m) 


JV-1 


1 

q(l,m)  =  —  ^  exp  (— i 2-Knl/N)p(n,m),  foil  =  Q:N  —  1,  m  =  0:M  —  1. 


n=0 


Therefore,  equation  (81)  for  the  average  error  becomes 

e  =  5>.("»)£|fl(l,*»)-ai(0  b(.l,m)\2. 


(82) 


(83) 


(84) 


To  minimize  e,  consider  the  partial  derivative 

-~YlWx (m) iq(k’ m) ~  ai 

1'  '  m 

for  k  =  ka  \ kb.  The  conditionally  optimal  coefficients  are  explicitly 


fli(*)  = 


q(k,m)  b*{k,m) 

Y^x(m)  |6(fc, m)|2 


for  k  =  ka.kb . 


(85) 


(86) 


These  solutions  are  decoupled  in  k  because  of  the  fiat  temporal  weighting  and  1/2"  frequency  separation. 
Substituting  equation  (86)  into  equation  (84)  yields  the  conditionally  minimal  error 

e  =  Yw*(m)J2 \-q(k’ m)  " fii w  m)3 9* (fc> m) 

m  Jfc 

=  E^(m)Ei5(fc>m)i2  -Eai(fc)Eu,*(m) 


=  YWx  E  i9(fc>  m)i2  -  E 


^^(m)  q{k,m)  b*(k,m) 


Ywx(.m)\b(k,m)\2 


(87) 
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Recall  (from  equations  (75)  and  (79))  that 
b(k,  m)  =  b(k,  m;  tq ,  1x2,  G,  A) 

=  ex.-p[-ia.kf}(rn)ui\  +  Gexpf-iajtA  -  iak0(m)u2]. 


(88) 


The  leading  term  in  equation  (87)  is  independent  of  hypothesized  values  iq,  112,  G,  and  A.  To  further 
minimize  e,  the  quantity 


r(ui,u2,G,X) 


M- 1 

^  wx(m)  q{k,m)  b*(k,m) 

m= 0 _ 

M  — 1 

53  w*(m)  l6(k>™)|2 

m=0 


(89) 


must  be  maximized  by  choice  of  ui,  U2,  G,  and  A.  The  sum  over  k  (frequency  bin  number)  reflects  an 
incoherent  addition,  because  no  phase  relationships  are  assumed  between  individual  frequency  components 
in  the  received  data.  The  numerator  sum  over  m  (spatial  element  number)  is  a  coherent  addition  when 
the  latter  four  parameters  of  b(k>m-,ui>U2,G,X)  line  up  with  the  actual  parameter  values  in  FFT  data 
{q(k,  m)}.  The  denominator  sum  over  m  in  equation  (89)  is  a  normalization  factor,  which  depends  on  k, 
the  frequency  bin  number,  in  addition  to  u\,  U2,  G,  and  A. 


Letting 

F(k)  =  y^wx ( m )  \q(k,  m) |2  for  k  =  ka:kb, 

m 

allows  equation  (87)  to  be  written  as 

4= (i-Pk), 


where 


Pk  = 


^2,wx(m)  b*(k,m)  q(k,m ) 


53 w*(m)  \Kk,m)\2^2wx(m)  |?(fc,m)|2 


By  Schwartz’s  inequality,  0  </?*,<  1  for  all  k. 


(90) 

(91) 


(92) 


Form  (92)  is  an  interesting  one,  but  equation  (89)  is  the  preferred  form  to  maximize.  The  “probe” 
function  b(k,m\ui,U2,G,  A)  in  equation  (88)  is  a  function  of  the  four  real  parameters  u  1,  U2,  G,  and  A; 
therefore,  a  four-dimensional  search  for  the  maximum  of  r  in  equation  (89)  is  required.  More  generally,  G 
could  be  complex,  if  desired;  then  a  five-dimensional  real  search  on  r  in  equation  (89)  would  be  needed. 


Once  the  optimal  parameter  values  (tti,  U2,  G,  A)  that  maximize  r  in  equation  (89)  have  been  deter¬ 
mined,  the  optimum  amplitudes  follow  from  equations  (85)  and  (78)  as 


ai (k)  = 


53wx(m)  q(k,m)  S*(A:,m) 

m 

53w*(m)  |f>(fc,m)|2 

771 


for  k  =  ka  '.kb, 


(93) 


where 


b(k,m)  =  exp[-iak/3(m)v,i]  +  Gexp[-zajfeA  -  iajfc/J(m)u2] .  (94) 

The  minimal  residual  in  the  frequency-space  domain  can  be  obtained  from  equation  (84)  in  the  form 

q(k,m)  =  q(k,m )  —  di(fc)  b(k,m).  (95) 
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5.  BEAMFORMING  FOR  A  SINGLE  MOVING  SOURCE 
NEAR  AN  ARBITRARY  PLANAR  ARRAY 

5.1  THEORY 


It  is  assumed  that  a  source  moves  from  (21,2/1)  to  (22,1/2)  during  an  observation  time  T,  with  constant 
speed  in  a  straight  line  (see  figure  22).  The  range  from  the  source  to  the  receiving  array  is  arbitrary;  also, 


(21,2/1  at  t  s=  0) 


_Source  motion 

*  (22,  at  t  =  T  —  NA) 


,y(m)}  for  m  =  0:M— 1 
- x 

Figure  22.  Sensor  and  Target  Motion 


the  range  and  angle  of  the  source  can  vary  during  the  observation  time  T.  Samples  are  taken  at  times 
nA,  where  n  =  0:  N  —  1  and  A  is  the  time-sampling  increment.  The  source  location  at  time  t  =  nA  is 


/  \  mL,2  —  •*'1 
xs(n)  =  xi+n  , 

,  2/2 -2/i 

ys{n)  =  2/1  +  n  N  _  1  • 


(96) 

(97) 


It  is  also  assumed  that  there  are  M  receivers  placed  at  coordinates  {x(m),y(m)},  with  m  —  0:M-1.  The 
separation  distance  between  the  single  source  and  the  m-th  receiving  element  at  time  nA  is 

s(n,m)  =  |[xs(n)  -  2(m)]2  +  [^(n) -y(m)]2|  .  (98) 

The  separation  distance  s(n,  m)  is  a  function  of  21,2/1,22,  2/2,  as  well  as  a  function  of  n  and  m. 

Modeling  the  source  as  transmitting  a  waveform  composed  of  a  sum  of  single  frequencies  yields  the 
corresponding  received  waveform  as: 

fcb 

53  a(k)  exp  (i  2-x  fkt),  (99) 

k=ka 


where  the  {a(k)}  are  complex  strength  coefficients.  The  fitting  frequencies  selected  will  be 


k 

NA 


for  k  =  ka:kb, 


(100) 
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although  no  simplification  results  from  this  choice.  Thus,  the  received  waveform  samples  are  modeled  at 
time  nA  and  element  m  as 


pi  (nA,i(m)) 


Y  a(k)  exp|i  2n  fk  [nA  —  ^s(n,  m)]  j 
k=ka 

Y  aC^expliafcfn-hCn.m)]},  a*  = 

k=zka  1 


(101) 


where  the  dimensionless  quantity 

h(n,  m)  =  ~~£~>  (102) 

and  c  is  the  propagation  speed. 


The  total  weighted  squared  error  between  this  model  and  the  actual  received  data  {p(n,  m)}  is 


JV-l  JW-l 


■-EE  wt(n)  wx (m) 


n= 0  m= 0 


kb 

' p(n,  m)  -  a(k)  exp  ji  ak  [n  -  h(n,  m)]  j 

k=k 


(103) 


where  a  set  of  (positive  real)  temporal  weights  {wt(n)}  and  (positive  real)  spatial  weights  {wx(m)}  are 
allowed.  For  fixed  hypothesized  x\,yi,X2,  and  y?,  the  set  of  amplitude  coefficients  {a(A:)}  that  minimizes 
error  e  is  to  be  found.  Finding  these  coefficients  requires  differentiation  of  error  e  with  respect  to  each  of 
the  amplitude  strengths,  indexed  by  k : 


de 


da*(k) 


where 


=  -  Y,wt(n)  ti>*(m)  p(n,  m)  —  Y,  °(fc)  expji  ak  \n  —  h(n,  m)]  j 

71)771  _  fc 

x  exp  j-i  a^[n  -  h(n,  m)]  j 

=  ~YWt  (n)  Wx  (m)  P(n> m)  exP  {  [n  ~  m)]  } 

n,m  '  * 

+  Y,  a(k)  Y,  wt  ( n )  wx  (m)  exp  |  —  (j:  —  A:^n  —  h(n,  m)  | 


k  n,m 

kb 


P(k) 
W(k ) 


=  -P(fc)  +  Y,  a(*0  W(k  —  k)  for  k  =  ka’-kb, 
k=ka 


Ywt(n )  wx(m)  p(n,m)  ex p  -  h{n,m )  j  , 

n,m  '  ' 

Y,  wtin)  w*(m)  exp  | ~i  [n  ~  m)]  |  • 

n  m  '  J 


(104) 


(105) 

(106) 


n,m 


The  last  quantity  is  defined  only  for  |/c|  <  kt,  -  ka.  W{k)  possesses  no  (5-function  properties,  regardless  of 
the  choices  of  weights  {u>t(n)}  and  {wx(m)}t  due  to  the  nonlinear  function  h(n,m). 
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The  conditionally  optimal  coefficients  {a (A:)}  are  found  by  setting  all  the  derivatives  in  equation  (104) 
equal  to  zero: 

kt 

53  W(k  -  k )  o(jfc)  =  P(4)  for  k  =  ka:kb.  (107) 

k=ka 

Equation  (107)  comprises  K  simultaneous  linear  equations,  where  K  =  kb  —  ka  +  1.  The  matrix  [W  (A:— A;)] 
is  Hermitian  Toeplitz  and  is  of  size  K  x  K.  The  corresponding  conditionally  minimal  error  e  is  calculated 
from  equations  (103)  and  (105): 


n,m  L  k 


=  Y^wt(n)  wx{m)\p(n,m)f  -  ^a(k)  P*(k).  (108) 

n,m  fc 


To  further  minimize  e  requires  maximization  of 

kt 

r(*i>  l/i »*2>  to)  =  53  flW  •?*(*)•  (109) 

k=ka 

Letting 

W  =[W(fc-Ac)],  a  =[fl(A:fl)..-a(i:t)]T,  P  =  [P(fca)  •  •  ■  P(kb)f ,  (110) 

A  XxV  Kx  1  A  XI 

results  in  equations  (107)  and  (109)  yielding 

W  a  =  P,  a  =  W_1  P,  r  =  P;  a  =  P'  W_1  P,  (111) 


where  the  prime  denotes  conjugate  transpose.  Vector  P  and  square  matrix  W  are  functions  of  x\t  j/i,  X2, 
and  2/2-  The  required  procedure  is  now  to  search  for  a  maximum  of  equations  (109)  or  (111)  in  Xi,y\,X2, 
and  2/2  space,  where  (g(A)}  are  the  solutions  to  equation  (107).  The  Hermitian  Toeplitz  character  of  W 
can  be  used  to  facilitate  solving  equation  (107),  rather  than  explicitly  calculating  its  inverse. 


Explicitly,  for  n  =  0:iV-l  and  m  =  0:M-1, 

1  2 


h(n,m)  = 


Xi+n^-~-x(m) 


,  2/2  -  2/i  ,  s 

m  +n  N~  i  -y(m> 


(112) 


For  each  hypothesized  x\,  y\,  x%,  and  z/2,  calculation  of  equation  (112)  requires  evaluation  of  NM  terms. 
Next,  equations  (105)  and  (106)  each  require  K  summations  of  size  NM.  Then,  equation  (107)  requires 
the  solution  of  K  simultaneous  linear  equations  with  a  Hermitian  Toeplitz  kernel.  Finally,  equation  (109) 
is  a  sum  of  K  complex  terms,  resulting  in  a  real  quantity. 


The  quantity  r(sBi,j/1,®2»I/2)  in  equation  (109)  is  maximized  at  position  estimates  ®i,  yi,  £2,  and  #2. 
First,  these  best  position  estimates  are  substituted  into  equation  (110),  and  matrix  W  and  vector  P  are 
evaluated.  Then,  Wa  =  P  is  solved  for  optimum  amplitudes  {o(A:)}.  Finally,  the  minimal  time-space 
residual  can  be  deduced  from  equation  (108)  in  the  form 

p(n,  m)  =  p(n,  m)  -  5 3  3(fc)  exp|i  a*  [n  -  h(n,  m)]  | ,  (113) 

where  h(n,  m)  follows  obviously  from  equation  (112).  There  is  no  simple  relation  for  the  minimal  frequency- 
space  residual  q(k,m )  because  of  the  nonlinear  dependence  of  h(n,m )  on  n. 
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5.2  IMPLEMENTATION 


The  theory  outlined  in  the  last  section  was  implemented  in  a  matlab  program  and  tested  on  simulated 
data  for  a  variety  of  scenarios. 


5.2.1  Data  Simulation 

The  first  step  was  to  generate  the  data  that  an  array  of  sensors  would  produce  in  the  presence  of  a 
moving  source.  The  starting  positions  and  ending  positions  were  defined  as  xl,  yi,  x2,  and  y2. 

y.starting  position: 
xl  =  110; 

yl  =  110; 

pending  position: 
x2  =  100; 

y2  =  120; 

The  basic  parameters  were  set  as  follows: 

c  =  1500;  '/.  wavespeed  (m/s) 

T  =20;  '/.  total  duration  of  data  (s) 

speed  =  5;  '/,  speed  of  source  (m/s) 

dt  =  .1;  '/.  sampling  increment  (s) 

ka  =  100;  '/.  starting  frequency  index 

kb  =  110;  '/,  ending  frequency  index 

M  =  8;  '/,  total  number  of  sensors 

dx  =  3.75;  '/,  spacing  of  hydrophones  (m) 

xarray  =  dx*(0:M-l) ; 
yarray  =  zeros (size (xarray) ) ; 

N  =  T/dt;  '/,  total  number  of  time  samples 
a  =  ones (size (ka: kb) ) ;  %  Frequency  components: 
d  =  speed*T;  ’/,  distance  traveled  by  source  (m) 

Some  variation  in  the  y  positions  of  the  receiving  array  elements  was  added  to  break  the  linearity  of  the 
array  and  the  resulting  left-right  ambiguity.  The  array  elements  were  shifted  in  a  random  way.  The 
random  number  generator  was  reset  to  its  initial  value  so  that  the  same  array  would  be  used  each  time 
the  program  was  run. 

randn(’ state' ,0) 
sigmay  =  1  '/,  (m) 

yarray  =  yarray  +  sigmay*randn(size(yarray)) ; 

Various  matrices  are  now  set  up  in  which  the  columns  correspond  to  the  different  elements  of  the  receiving 
array,  and  the  rows  correspond  to  the  different  time  samples.  First,  “matrix”  versions  of  the  element 
locations  are  set  up,  as  follows: 
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xm  =  ones(N,l)*xarray; 
yns  =  ones(N, l)*yarray ; 


Then,  matrices  of  indices,  which  will  be  used  later,  are  set  up.  These  matrices,  m  and  n,  are  arranged  such 
that  the  elements  of  m  vary  only  across  its  columns  and  the  elements  of  n  vary  only  down  its  rows: 


Qn,n]  =  meshgrid(0: (M  -  1),0:(N  -  1)); 


For  example,  the  first  few  elements  of  m  and  n  are 


»  m(l:3,l:3) 
ans  = 

0  12 

0  12 

0  12 

»  n(l:3,l:3) 
ans  = 

0  0  0 

111 
2  2  2 


These  matrices  can  then  be  used  in  a  “vectorized”  sense  to  calculate  quantities  such  as  the  source  location 
at  each  timestep: 


xsn  =  xl  +  n*(x2  -  xl)/(N  -  1) ; 
ysn  =  yl  +  n*(y2  -  yl)/(N  -  1) ; 


The  dimensionless  separation  distance  between  the  source  and  the  m-th  receiving  element  at  timestep  n 
can  be  calculated  according  to: 


hum  =  sqrt((xsn  -  xm).~2  +  (ysn  -  ym) .  "2)/(c*dt) ; 


This  calculation  corresponds  to  equation  (112).  A  f  or-loop  is  used  to  add  the  frequency  components  to 
obtain  the  signal  received  by  each  element  at  each  time  step: 

pnm  =  0; 
k  =  ka:kb; 
i  *  sqrt(-l) ; 
twopioverN  =  2*pi/N; 
for  this_k  =  1: length (k) 

alphak  =  twopioverN*k(this_k) ; 
pnm  =  pnm  +  a (this _k)*exp (alphak* (n  -  hnm)*i); 
end 
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The  matrix  pnm  is  the  data  matrix,  in  which  time  runs  down  the  rows  and  the  hydrophone  number  runs 
across  the  columns. 

Next,  the  option  of  adding  some  complex  normally-distributed  random  noise  using  the  randn  function 
is  available: 


snr  =  0;  V,  dB 

Noise  =  randn  (size  (pnm))  +  i*randn  (size  (pnm) ) ; 
As  =  mean (std (pnm)) ; 

An  =  mean(std(Noise)) ; 

NoiseScaleFactor  =  As  /  (An*10"(snr/20)) ; 
pnm  =  pnm  +  Noise  *  NoiseScaleFactor; 


5.2.2  Example 


An  example  of  the  output  produced  by  the  simulation  coding  is  shown  in  figure  23,  where  a  linear 
receiving  array  of  M  =  10  sensors  is  seen  in  the  top  plot.  The  source  emitted  a  single  frequency  tone,  and 
moved  from  the  beginning  of  the  arrow'  to  the  end  of  the  arrow  during  the  acquisition  time.  The  lower 
two  plots  show  the  signal  received  at  each  of  the  elements  during  a  short  period  at  the  beginning  of  the 
acquisition  time  (left)  and  at  the  end  (right).  The  light  and  dark  strips  represent  the  peaks  and  troughs 
of  the  received  waveform.  At  the  beginning  of  the  acquisition  time,  the  curved  wavefronts  arrive  at  the 
elements  near  the  middle  of  the  array  first.  The  wavefronts  can  be  imagined  as  a  series  of  concentric  circles 
emanating  from  the  tail  of  the  arrow.  At  the  end  of  the  acquisition  time,  the  curved  wavefronts  arrive 
at  element  10  first  and  then  propagate  down  the  array  towards  element  1.  In  this  case,  the  wavefronts 
can  be  imagined  as  a  series  of  concentric  circles  emanating  from  the  head  of  the  arrow.  An  animation 
of  the  strip  display  would  show,  at  the  start  of  the  acquisition,  the  wavefronts  meeting  the  array  first  at 
the  middle  of  the  array  (plot  at  left).  Then,  as  time  elapsed,  the  wavefronts  would  meet  the  array  at 


y 

source 


[1  23456789  10 


receiving  array 
—  x 


0  0.02  0.04 

Time,  s 


4.95  4.96  4.97  4.98  4.99 
Time,  s 


Figure  23.  Simulated  Data 
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points  progressively  closer  to  element  10  (plot  at  right).  The  optimization  algorithm  amounts  to  fitting  a 
modeled  sum  of  frequencies  to  these  data  by  choosing  the  best  start  and  end  points  for  the  arrow  (xi,  y i, 
x2,  and  j/2). 


5.3  OPTIMIZATION  ALGORITHM 


The  driver  code  for  the  optimization  used  the  matlab  function  fminsearch  to  do  the  minimization. 
This  function  employed  the  Nelder-Mead  simplex  (direct  search)  method  for  multidimensional,  uncon¬ 
strained  nonlinear  minimization.7  The  algorithm  returns  a  vector  that  is  a  local  minimizer  of  the  function 
movingsourceerror  near  the  starting  vector  guess: 


guess  =  [xguessl  yguessl  xguess2  yguess2] ; 

opts  =  optimset( ’Maxlter ’  ,2000,  ’MaxFunEvals ’  ,2000) ; 

x  =  fminsearch( ’movingsourceerror ’  .guess, opts) ; 


The  variable  opts  controls  the  maximum  allowable  number  of  iterations  and  function  evaluations. 

The  function  movingsourceerror  contains  the  code  that  first  finds  the  optimal  coefficients  ak  by 
solving  the  simultaneous  equation  (107),  and  then  uses  these  coefficients  to  calculate  the  quantity 
f(xi,yi,X2,y2)  given  by  equation  (109).  The  code  begins  by  defining  the  time  and  space  weighting 
functions  as: 


wt  =  ones(l,N)/N; 
wx  =  ones (1 ,M)/M; 

It  then  computes  the  quantity  h(n,m )  based  on  the  parameters  Xj,  j/i,  x2,  and  i/2,  which  are  given  as 
input  arguments  to  the  function: 


hnm  =  sqrt((xl  +  n*(x2  -  xl)/(N  -  1)  -  xm)  ."2  +  ... 

(yl  +  n*(y2  -  yl)/(N  -  1)  -  ym) . "2) /(c*dt) ; 


The  matrix  W  is  a  function  of  (k-k),  with  a  sum  over  m  and  n  required  to  calculate  each  element. 
Because  all  the  dimensions  for  standard  Matlab  two-dimensional  matrices  have  been  used,  one  must 
index  over  n  and  m  and  use  a  matrix  to  represent  (k  —  k).  Two  f  or-loops  can  then  be  used  to  perform 
the  summation  over  n  and  m: 


kbar  =  k; 

[K.Kbar]  =  meshgrid(k,kbar) ; 
W  =  0; 

P  =  0; 

KbarminusK  =  Kbar  -  K; 


Since  there  is  no  zero  index  in  Matlab,  loops  going  from  nn  =  1:N,  etc.,  are  used  as  matrix  indices, 
whereas  nn-1  etc.,  are  used  in  the  calculations: 
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for  mm  =  1:M 
for  nn  =  1:N 

fact  =  -2*pi*((nn  -  1)  -  hnm(im,mm))/N; 
weight  =  wt (nn) *wx (mm) ; 

W  =  W  +  weight*exp(fact*KbarminusK*i) ; 

P  =  P  +  weight*exp(fact*k*i)*pnm(nn,mm) ; 
end 
end 


The  theory  above  has  P  as  a  column  vector;  therefore,  the  transpose 
P  =  P.»; 

is  chosen.  Now  the  equation  W A  =  P  is  solved  to  obtain  the  least-squares  coefficients,  {a(fc)},  which  will 
be  returned  in  the  vector  A: 

A  =  W\P; 

Having  found  the  conditionally  optimal  {a (A:)},  r  must  be  maximized  to  minimize  the  error  e.  However, 
because  a  minimizer  is  being  used,  the  sign  must  be  changed: 

r  =  -real(P’*A) ; 

(The  imaginary  part  of  r  in  equation  (109)  should  be  zero;  therefore,  the  smallness  of  the  magnitude  of 
the  imaginary  part  of  r  forms  a  partial  check  on  the  accuracy  of  the  final  solution.) 


5.4  RESULTS 

5.4-1  Examination  of  the  Error  Function 

The  minimal  error  e  given  in  equation  (108)  can  be  considered  as  a  scalar  function  of  (xi, yi, £2,2/2)- 
To  understand  the  behavior  of  this  (four-dimensional)  function,  e  is  calculated  as  a  function  of  either  the 
start-point  coordinates  with  the  end  point  fixed,  or  as  a  function  of  the  end-point  coordinates  with  the 
start  point  fixed,  both  of  which  are  two-dimensional  “slices”  through  the  four-dimensional  function.  Two 
cases  are  considered:  in  one,  the  source  moves  in  the  nearfield  of  the  array,  and  in  the  other,  the  source 
moves  in  the  farfield  of  the  array. 
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5.4-1. 1  Nearfield  Scenario.  Figure  24  shows  the  variation  in  the  conditionally  minimal  error  e 
with  changes  in  the  modeled  start  and  end  points  of  the  trajectory.  The  upper  plots  show  contours  of 
the  error  quantity  e.  The  lower  plots  show  the  same  quantity  plotted  as  surfaces.  Two  scenarios  are 
illustrated.  The  plots  on  the  left  (top  and  bottom)  show  e  as  a  function  of  (21,2/1)  with  X2  and  2/2  held 
fixed  at  their  true  values,  and  the  plots  on  the  right  (top  and  bottom)  show  e  as  a  function  of  (22,2/2)  with 
X\  and  j/i  held  fixed  at  their  true  values.  Contours  are  drawn  at  the  same  levels  in  both  contour  plots. 
The  dots  show  the  positions  of  the  sensor  array,  and  the  arrow  shows  the  true  trajectory  of  the  source 
during  the  observation  period.  In  both  scenarios,  there  is  a  deep  minimum  in  the  region  within  roughly 
a  wavelength  of  the  true  value.  (These  runs  simulated  a  noise-free,  single  planewave,  with  k  =  400,  and 
T  =  5  s,  yielding  a  wavelength  of  18.75  m.)  Although  a  line  array  receiver  was  used,  the  resulting  e  is  not 
expected  to  be  left-right  symmetric  because  the  other  end  of  the  trajectory  was  held  fixed;  therefore,  the 
left-right  symmetric  point  would  represent  a  completely  different  source  trajectory  that  is  easily  ‘rejected’ 
by  the  fitting  process.  The  error  function  would  be  left-right  symmetric  if  either  end  of  the  trajectory  lay 
along  the  axis  of  the  array.  Outside  the  region  within  a  wavelength  of  the  true  solution,  the  error  function 
is  relatively  flat  and  appears  to  contain  many  shallow  local  minima.  This  kind  of  function  may  present  a 
challenge  to  the  minimization  algorithm  chosen,  especially  in  the  presence  of  noise. 


Figure  24  ■  Error  for  Nearfield  Scenario 


5.4-1-2  Farfield  Scenario.  Figure  25  shows  the  same  type  of  results  as  figure  24,  except  this  time 
the  source  is  at  a  greater  distance  from  the  array.  The  error  function  has  a  narrow  valley  that  follows  an 
arc  centered  near  the  location  of  the  array.  Again,  this  kind  of  function  may  present  a  challenge  to  the 
minimization  algorithm  chosen,  especially  in  the  presence  of  noise. 
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Figure  25.  Error  for  Farfield  Scenario 


5-4.2  Fitted  Trajectories  with  Noisy  Data 

Simulations  in  this  section  were  carried  out  using  the  scenario  shown  in  figure  26.  An  array  of  eight 
elements  w*as  located  along  the  z-axis,  with  each  element  equally  spaced  at  a  separation  distance  of 
3.75  m.  The  first  element  of  the  array  was  located  at  the  origin.  A  source  moving  in  the  nearfield  of 
the  array  (as  shown  by  the  thick  arrow  indicated  on  the  plot)  emitted  frequencies  /*  =  k/T,  where 
k  =  ka'-k/,,  ka  =  100,  kb  =  110,  and  T  —  20  s.  The  speed  of  the  source  was  such  that  it  moved  from 
(xi,2/i)  =  (60,60)  to  (22,3/2)  =  (50,70)  during  the  time  interval  T  =  20  s.  The  propagation  speed  was 
taken  to  be  1500  m/s.  The  time  between  samples  was  A  =  0.1  s,  resulting  in  N  =  200  samples  per 
element.  Noise  was  added  to  the  pressure  data  at  various  SNRs;  the  SNR  is  defined  in  the  same  way  as  on 
page  42.  A  search  was  carried  out  for  the  values  of  21, 3/1 , 22,  and  2/2  that  minimized  e,  as  described  above. 
Five  realizations  of  the  random  noise  field  were  used  at  each  SNR.  The  values  used  to  initialize  the  search, 
which  remained  constant  throughout  these  runs,  consisted  of  motion  beginning  at  a  range  of  500  m  at  a 
bearing  equal  to  tan-1  { [(3/1  +  3/2 ) /2]  /  [(zi  +  Z2)/2]  }  and  finishing  at  the  position  reached  after  traveling 
for  time  T  at  10  knots  towards  the  origin.  The  numerical  values  were  xi  =  323.0, 3/1  =  382.0,  x2  =  257.0, 
and  3/2  =  303.0.  The  results  of  the  minimizations  are  shown  by  the  thin  arrows,  with  the  size  of  the 
arrow-heads  being  inversely  proportional  to  the  SNR,  as  indicated  on  the  legend. 

The  data  of  figure  26  are  presented  as  individual  coordinate  errors  and  as  functions  of  SNR  in  figure  27. 
The  four  plots  show  the  magnitude  of  the  difference,  in  meters,  between  the  coordinates  obtained  from 
the  search  and  the  true  coordinates,  as  a  function  of  SNR.  The  results  show  that,  even  for  this  close-range 
case  (85  m),  an  error  on  the  order  of  10  m  remains  for  SNRs  less  than  about  40  dB. 


46 


SNR,  dB 
&►() 

►  10 
e-  20 
■®*-  30 

—  40 
-*►50 

—  60 

—  70 

—  80 

—  90 

—  100 


Figure  26.  True  Trajectory  and  Fitted  Trajectories  for  Scenarios  with  Various  SNRs 
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Figure  27.  Fitted  Coordinate  Errors  Plotted  Against  SNR 
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5.4-3  Comparison  with  Conventional  Beamformer 


The  new  algorithm  has  the  potential  to  outperform  beamformers  because  it  yields  tracking  information 
directly:  in  conventional  sonar  systems,  tracking  is  performed  only  after  a  set  of  peaks  in  the  beamformed 
data  has  been  obtained.  But,  for  a  moving  source,  such  peaks  may  never  be  detected  because  the  beam- 
former  “peak”  will  be  smeared  out  over  the  range  of  angles  swept  out  by  the  moving  source  during  the 
integration  time  of  the  beamformer.  An  example  of  this  can  be  seen  in  figure  28.  The  plots  at  the  top 
are  for  a  stationary  source,  and  the  plots  at  the  bottom  are  for  a  moving  source.  The  beamformed  output 
for  the  stationary  source  shows  a  single  peak  at  the  correct  bearing,  but  the  beamformed  output  for  the 
moving  source  shows  energy  spread  out  over  the  angles  swept  out  by  the  source  as  it  moves  during  the 
time  over  which  the  data  were  acquired  (depicted  by  the  gray  areas).  The  output  in  this  case  could  be 
mistaken  for  a  pair  of  stationary  sources.  It  is  also  possible  that  an  automatic  detector  would  not  yield  a 
detection  for  the  moving  case  because  of  its  high  sidelobe  level. 


Source^ 


Stationary1 

source 


20  metres 

i - 1 


Moving  source 


Conventional  beamformer  results 


Figure  28. 


Conventionally  Beamformed  Results  for  a  Stationary  Source  (Top) 
and  a  Moving  Source  (Bottom) 
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6.  CONCLUSIONS 


A  new  array  processing  technique  has  been  presented  that  allows  optimal  source  amplitudes  and  arrival 
angles  to  be  determined  from  data  produced  by  arbitrary  arrays.  Used  to  fit  either  one  or  two  planewaves 
to  arbitrarily  spaced  linear  array  data,  the  technique  has  shown  superior  performance  to  conventional 
beamforming  when  more  than  one  source  is  within  the  beamwidth  of  an  array.  With  such  situations 
becoming  more  likely  as  operating  frequencies  are  driven  to  lower  values,  in  an  effort  to  increase  detection 
ranges,  the  approach  can  be  used  to  coherently  subtract  strong  arrivals  from  received  data,  allowing  the 
detection  of  additional  weak  arrivals  in  close  angular  proximity. 

When  a  planewave  arrives  at  an  arbitrary  two-dimensional  array  from  a  moving  source,  a  similar 
technique  yields  optimal  values  of  the  source  amplitudes  of  the  planewave,  as  well  as  the  starting  and 
ending  positions  of  the  moving  source.  The  new  processing  thus  performs  the  functions  of  a  combined 
beamformer  and  tracker  in  such  a  way  that  the  total  error  over  the  array  and  the  observation  interval  is 
minimized. 

The  results  presented  here  are  for  simulated  data.  Future  work  will  involve  the  application  of  these 
methods  to  measured  at-sea  data.  A  further  issue  to  be  addressed  is  the  adequacy  of  the  four-dimensional 
iterative  minimization  procedures. 
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APPENDIX 


ALTERNATIVE  FORMULATION 
OF  TWO-PLANEWAVE  FIT 

There  is  an  alternative  formulation  to  the  two-planewave  fitting  waveform  in  equation  (46)  that  is  exact 
for  all  U\,U2  and  does  not  develop  any  singularities  as  u\  -¥  u2.  Letting  central  angle  uc  and  difference 
angle  e  be  defined  as 


uc  =  (ui+u2)/2,  e  =  ui  -  u2, 


(A-l) 


in  terms  of  the  original  arrival  angles  Hi  and  u2  (difference  angle  e  need  not  be  small),  then,  by  using 
u1=uc  +  e/2,  u2  =  uc  —  e/2,  the  k- th  frequency  component  of  the  original  model  (equation  (46))  can  be 
expressed  as 


th  =  bi(k)  Bi(k,m)  +  b2(k)  B2(k,m)  exp (iakn)  for  k  =  ka  -.kb, 


(A-2) 


where 


2irk  „  x(m)  .  .  .  sin(rr) 

“*  =  -jf,  0m  =  ,  sm*(x)  =  — , 

Bi(k,m )  =  cos(a*/?me/2)  exp(-i  akpmuc) , 
B2(k,m )  =  sinx(a*:/3me/2)  exp(-iaA/3muc), 

b\(k)  =  a\(k)  +  a2(k), 
b2(k)  =  -i  [aj(A:)  —  a2{k)\ake/ 2. 


(A-3) 


For  the  best  fit  of  sum  waveform  Y,k  **  t0  a  given  data  set  {p{n,  m)},  the  optimal  coefficients  {i>i  (A:)}  and 
{b2(k)}  in  equation  (A-2)  will  all  be  finite,  regardless  of  the  size  of  e,  including  e  =  0.  However,  because 


ai(k)  =  §&i(A;)  + 


ib2(k) 


®a(*)  =  %h(k)  - 


ib2(k) 
ake  ’ 


(A-4) 


then  di(k)  — >  oo  and  a2{k)  — y  oc  as  e  0  in  the  original  model  (equation  (46)).  (This  effect  has 
been  observed  before  and  has  been  explained  by  Nuttall.8)  On  the  other  hand,  the  model  formulation  of 
equation  (A-2)  does  not  undergo  this  singular  behavior  for  any  e. 


As  e  -»  0,  the  model  waveform  tk  becomes 
^2  [6i  (k)  +  b2  (k)  Pm]  exp  [i  ak  (n  -  Pmuc)] , 


(A-5) 


upon  use  of  equations  (A-2)  and  (A-3).  This  dependence  of  the  fit  on  m  is  different  from  the  original 
model  (equation  (46));  such  modified  behavior  has  been  described  in  eq.  28  et  seq.  and  eq.  51  et  seq.  in 
Nuttall.8 

When  component  form  (A-2)  is  used  to  approximate  the  available  data,  the  error  to  be  minimized 
becomes,  for  flat  temporal  weighting, 


=  Iy 

N 


kb 


\p(n,m)  -  *k 


k—ka 


(A-6) 


A-l 


Holding  angles  uc  and  e  fixed  for  now,  the  partial  derivatives  of  e  with  respect  to  {bi(k)}  and  {bi(k)} 
yield  decoupled  pairs  of  simultaneous  linear  equations  for  the  conditionally  best  coefficients  {bi(k)}  and 
{b2(k)}  as 

WnWb^  +  WuWb^^Qiik)  ) 

>  for  k  =  ka:kb,  (A- 7) 

W21(k)  b, (k)  +  W22(k)b2(k)  =  Q2(k)  J 


where 


Wn(k)  =  'Y^Wxim)  |Bi(A:,m)|2  =  £^3(77-1)  cos2 (ak/3mef 2), 

m  TTi 

W22(k)  -^2wx(m)\B2(k,m)\2  =  ^ w* (m)^ sinx2(a*/3m<:/2), 

m  m 


Wi2(k)  =  W21  (k)  =  ^2wx(m)B*1(k,m)B2(k,Tn) 

m 

=  £  wx(m)0m  sinx(afc/3me/2)  cos(ajt/3me/2), 

m 

Qj{k)  =  £tw*(m)0(*>m)Bi(*:«m)>  j  =  1,2. 


m 


\ 


> 


j 


(A-8) 


It  should  be  observed  that  all  the  (Wjj(A:)}  are  real  functions  of  e,  but  are  independent  of  uc  and  the  data 
{p(n,m)}  or  {q(k,m)}.  The  quantities  {Bj(k,m)}  are  complex  functions  of  e  and  uc,  but  are  independent 
of  the  data.  The  quantities  {<2j(fc)}  are  complex  functions  of  e  and  uc  as  well  as  complex  functions  of  the 
data. 


The  solutions  to  equation  (A-7)  are  explicitly 


,  (1,  _W22{k)Qx{k) -W12(k)Q2{k)  \ 

-lU  Wn(k)W22(k)-W?2(h)  I 

)  for  k  —  fcft  i  kh 

.  _  W11(k)Q2(k)-W12(k)Q1(k) 

22  u  W11(k)W22(k)-W^(k)  J 


(A-9) 


In  the  limit  as  e  -»  0,  it  follows  from  equation  (A-8)  that  the  denominator  of  equation  (A-9)  tends  to 


<  m  /  \  m  J 


(A-10) 


which  is  always  positive  by  Schwartz’s  inequality.  Thus,  equation  (A-9)  never  develops  any  singularities 
for  any  value  of  e;  that  is,  the  denominator  of  equation  (A-9)  is  never  zero  for  any  e. 


When  the  conditionally  optimal  coefficients  in  equation  (A-9)  are  substituted  in  error  expression  (A-6), 
and  equations  (A-2)  and  (A-8)  are  employed,  the  conditionally  optimal  error  becomes 


e=  ^£w.r(m)  p(n,m)  —  £{Mfc)si(*;>7n)  +h(k)B2(k’m)} 


=  ^£«>x(m)b(n,m)|2  -  £{m*)QI(*)  +  b2(k)Q'2(k)} 

nro  k 

=  jj  £  Wx (m) |p(n, m) |2  -r2(uc,e), 

nm 


p*  (n,  m) 


(A-ll) 

(A-12) 

(A-13) 


A-2 


where  function 


r  ,  v  _  Wn{k)\Qi{k)?  +  Wu(k)\Q2(k)\2  -  2W12(W{Ql(k)Q2(h)} 

Wn(k)WxL(k)-W&(k) 


(A-14) 


and  K  denotes  the  real  part.  To  further  minimize  conditional  error  e,  the  quantity  r2(uc,e)  must  now 
be  maximized  by  choice  of  central  angle  uc  and  difference  angle  e;  see  equation  (A-l).  No  singularities 
develop  in  r2(uc,e)  for  any  value  of  e,  including  e  =  0. 


Consideration  of  equation  (A-8)  reveals  that  {W„(A:)}  are  real  functions  of  e  and  k;  therefore,  their 
values  could  be  precalculated  and  stored,  at  least  for  a  preliminary  coarse  search  where  difference  angle  e 
could  be  fixed.  Similarly,  the  quantities  {Bj(k,m)}  in  equation  (A-3)  are  complex  functions  of  e  and  uc, 
but  are  independent  of  the  data  (p(n,m)}  or  {q(k, m)}.  For  a  fixed  e,  {Bj(k,m)}  could  also  be  stored; 
the  size  of  the  storage  would  be  KMMC,  where  K  is  the  number  of  frequency  components  in  the  original 
model  (equation  (46)),  M  is  the  number  of  array  elements,  and  Mc  is  the  number  of  central  angles  uc  of 
interest.  On  the  other  hand,  the  quantities  {Qj(k)}  in  equation  (A-8)  are  complex  functions  of  e  and  uc, 
as  well  as  complex  functions  of  the  data  {q(k,  m)},  and  cannot  be  precalculated. 

Once  the  unconditionally  optimal  angles  uc  and  i  are  found  from  the  maximization  of  r2(uc,e)  in 
equation  (A-14),  these  angles  must  be  substituted  into  equation  (A-3)  to  obtain  {Bj(k,  m)},  and  then  into 
equation  (A-8)  to  obtain  {Wij(k)}  and  {(?j(^)}-  Finally,  equation  (A-9)  yields  the  best  coefficient  values 
{fey  (A) } ,  where  all  the  quantities  must  have  their  optimal  (*)  values. 

The  minimal  residual  in  the  time-space  domain  follows  from  equations  (46)  and  (A-2)  as 
p(n,m)  =  p(n,m) -p2(nA,x(m))  (A-15) 

=  p(n,  m)  -  ^2  J^fei  ( k)Bi  (k,  m,  e,  uc)  +  S2  (k)B2(k,  m,  e,  uc)j  exp(i  ajtn)  (A-16) 

for  n  —  0:N—1,  m  =  0:Af— 1.  The  corresponding  minimal  residual  in  the  frequency-space  domain  is 
simply 

q(k,m)  =  q(k,m)  -bi(k)Bi(k,m,i,uc)  -b2{k)B2(k,m,e.,uc)  (A-17) 

for  k  =  ka:kb,  m  =  0:M-1.  These  two  forms  of  coherent  subtraction  should  be  employed  using  coeffi¬ 
cients  {bj(k)}  rather  than  resorting  to  the  original  model  (equation  (46))  and  coefficients  {aj(k)}  from 
equation  (A-4).  The  latter  coefficients  will  be  extremely  large  for  very  small  e  and  will  cause  loss  of 
significance  in  residuals  {p(n, m)}  and  {q(k, m)}. 

If  e  is  not  very  small,  the  optimal  amplitudes  in  the  original  model  (equation  (46))  can  be  found  from 
equation  (A-4)  as 


di(fc)  =  hti(k)  +  a2{k)  =  iSi(fc)  -  iMQ.  (A.18) 

QfcC  Qk€ 

This  relation  is  exact  and  applies  for  all  e  ^  0. 

This  appendix  has  presented  an  alternative  approach  and  solution  to  the  two-planewave  fitting  proce¬ 
dure;  all  quantities  are  finite  for  all  e,  except  possibly  for  equation  (A- 18),  which  attempts  a  return  to  the 
amplitudes  in  the  original  model  (equation  (46)). 
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